%~Mouliné par MaN_auto v.0.37.1 (34d45bef) 2025-11-24 15:13:53
\documentclass[OJMO,Unicode,published]{cedram}

\usepackage[all]{xy}
\usepackage{csquotes}
\usepackage{acro}
\usepackage{nicefrac}
\usepackage[nameinlink]{cleveref}
\usepackage{mathtools}
\usepackage{tabularx}
\usepackage{booktabs}
\usepackage{mleftright}

\crefname{assu}{Assumption}{Assumptions}
\crefname{theo}{Theorem}{Theorems}
\crefname{lemm}{Lemma}{Lemmata}
\crefname{coro}{Corollary}{Corollaries}
\crefname{rema}{Remark}{Remarks}
\crefname{defi}{Definition}{Definitions}
\crefname{section}{Section}{Sections}
\crefname{subsection}{Section}{Sections}
\crefname{equation}{equation}{equations}
\crefname{figure}{Figure}{Figures}

\crefname{step}{step}{steps}

\newcommand{\CC}{\mathtt{c}}

\newcommand*{\set}[1]{\mleft\{#1\mright\}}
\newcommand*{\args}[1]{\mleft(#1\mright)}

\newcommand*{\compromise}{\texttt{Compromise.jl}}
\newcommand*{\testmo}{\texttt{TESTMO}}
\newcommand*{\dfmo}{\texttt{DFMO}}
\newcommand*{\cset}{\texttt{CSET}}
\newcommand*{\cseq}{\texttt{CSEQ}}
\newcommand*{\cpen}{\texttt{CPEN}}

\newcommand*{\solvers}{\mathbf{Solvers}}
\newcommand*{\problems}{\mathbf{Problems}}

\newcommand*{\tangcone}[2]{\mathcal T_{#1}\args{#2}}
\newcommand*{\normcone}[2]{\mathcal N_{#1}\args{#2}}

\newcommand{\abs}[1]{\left| #1 \right|}
\newcommand*{\norm}[1]{\left\|#1\right\|}

\newcommand*{\ve}[1]{\pmb{#1}}
\newcommand*{\mat}[1]{\ve{#1}}

\newcommand*{\eye}[1]{\mat{I}_{#1,#1}}

\newcommand{\grad}{\ve{\nabla}\!}
\newcommand{\jac}{\grad}
\newcommand{\hess}{\ve{\nabla}^2\!}

\newcommand*{\tran}[2][2]{#2^{\mkern-#1mu\mathsf{T}}}

\newcommand*{\dotProd}[3][2]{\tran[#1]{#2}\!\cdot\!#3}
\newcommand*{\noDotProd}[3][2]{\tran[#1]{#2}\,#3}
\newcommand*{\dotMatProd}[2]{#1\!\cdot\!#2}
\newcommand*{\noDotMatProd}[2]{#1\,#2}

\newcommand*{\kk}{k}

\newcommand*{\iterUpCmd}[1]{\ensuremath{^{(#1)}}}
\newcommand*{\iterUp}{\iterUpCmd{\kk}}
\newcommand*{\iterDownCmd}[1]{\ensuremath{_{(#1)}}}
\newcommand*{\iterDown}{\iterDownCmd{\kk}}
\newcommand*{\iterSubCmd}[1]{\ensuremath{_{#1}}}
\newcommand*{\iterSub}{\iterSubCmd{\kk}}

\newcommand*{\itCmd}[1]{\iterUpCmd{#1}}
\newcommand*{\itDeltaCmd}[1]{\iterDownCmd{#1}}
\newcommand*{\itthetaCmd}[1]{\iterSubCmd{#1}}

\newcommand*{\KK}{\itCmd{\kk}}
\newcommand*{\KKn}{\itCmd{\kk + 1}}

\newcommand*{\normTr}[1]{\norm{#1}_{\mathrm{tr}}}
\newcommand*{\normTrIt}[1]{\norm{#1}_{\mathrm{tr,\kk}}}
\newcommand*{\normTwo}[1]{\norm{#1}_2}
\newcommand*{\normInf}[1]{\norm{#1}_\infty}
\newcommand*{\normIt}[1]{\norm{#1}\iterSub}

\newcommand*{\feas}{\mathcal X}
\newcommand*{\vx}{\ve{x}}
\newcommand*{\vZ}{\ve{0}}

\newcommand*{\vd}{\ve{d}}
\newcommand*{\vs}{\ve{s}}
\newcommand*{\vt}{\ve{t}}
\newcommand*{\vn}{\ve{n}}
\newcommand*{\vv}{\ve{v}}
\newcommand*{\vr}{\ve{r}}
\newcommand*{\vdk}{\vd\KK}
\newcommand*{\vnk}{\vn\KK}
\newcommand*{\vnkN}{\vn\KKn}
\newcommand*{\vsk}{\vs\KK}
\newcommand*{\vtk}{\vt\KK}
\newcommand*{\vvk}{\vv\KK}
\newcommand*{\vrk}{\vr\KK}
\newcommand*{\vxNk}{\vxk_{\mathrm{n}}}
\newcommand*{\vxN}{\vx_{\mathrm{n}}}
\newcommand*{\vxSk}{\vxk_{\mathrm{s}}}
\newcommand*{\vxS}{\vx_{\mathrm{s}}}

\newcommand{\vxk}{\vx\KK}

\newcommand*{\fL}{f}
\newcommand*{\gL}{g}
\newcommand*{\hL}{h}

\newcommand*{\fLk}{\fL\KK}
\newcommand*{\gLk}{\gL\KK}
\newcommand*{\hLk}{\hL\KK}

\newcommand*{\surrogateSign}{\hat}

\newcommand*{\vf}{\ve{\fL}}
\newcommand*{\vh}{\ve{\hL}}
\newcommand*{\vg}{\ve{\gL}}

\newcommand*{\indf}{\mathfrak{l}}
\newcommand*{\indg}{\mathfrak{i}}
\newcommand*{\indh}{\mathfrak{j}}
\newcommand*{\ind}{\mathfrak{m}}

\newcommand*{\indmax}[1]{#1_{\!\text{\scriptsize*}}}

\newcommand*{\df}{\grad{\fL}}
\newcommand*{\dg}{\grad{\gL}}
\renewcommand*{\dh}{\grad{\hL}}

\newcommand*{\dfk}{\df\KK}
\newcommand*{\dgk}{\dg\KK}
\newcommand*{\dhk}{\dh\KK}

\newcommand*{\DF}{\mat{\MakeUppercase{\fL}}}
\renewcommand*{\DH}{\mat{\MakeUppercase{\hL}}}
\newcommand*{\DG}{\mat{\MakeUppercase{\gL}}}

\newcommand*{\DFk}{\DF\iterSub}

\newcommand*{\varDim}{n}
\newcommand*{\fDim}{K}
\newcommand*{\hDim}{M}
\newcommand*{\gDim}{P}

\newcommand*{\mfL}{\surrogateSign{\fL}}
\newcommand*{\mhL}{\surrogateSign{\hL}}
\newcommand*{\mgL}{\surrogateSign{\gL}}

\newcommand*{\mf}{\ve{\surrogateSign{\fL}}}
\newcommand*{\mg}{\ve{\surrogateSign{\gL}}}
\newcommand*{\mh}{\ve{\surrogateSign{\hL}}}

\newcommand*{\mfLk}{\mfL\KK}
\newcommand*{\mgLk}{\mgL\KK}
\newcommand*{\mhLk}{\mhL\KK}

\newcommand*{\mfk}{\mf\KK}
\newcommand*{\mgk}{\mg\KK}
\newcommand*{\mhk}{\mh\KK}

\newcommand*{\mgkSub}{\mg\iterSub}

\newcommand*{\dmf}{\grad{\mfL}}
\newcommand*{\dmg}{\grad{\mgL}}
\newcommand*{\dmh}{\grad{\mhL}}
\newcommand*{\dmfk}{\dmf\KK}
\newcommand*{\dmgk}{\dmg\KK}
\newcommand*{\dmhk}{\dmh\KK}


\newcommand*{\DmF}{\mat{\surrogateSign{\MakeUppercase{\fL}}}\iterSub}
\newcommand*{\DmH}{\mat{\surrogateSign{\MakeUppercase{\hL}}}\iterSub}
\newcommand*{\DmG}{\mat{\surrogateSign{\MakeUppercase{\gL}}}\iterSub}

\newcommand*{\scal}[1]{\Phi\bigl[#1 \bigr]\!}
\newcommand*{\scalf}[1]{\scal{\vf}\args{#1}}
\newcommand*{\scalShort}{\Phi}
\newcommand*{\scalk}{\scalShort\KK}

\newcommand*{\stepsize}{\sigma}
\newcommand*{\stepsizek}{\stepsize\iterDown}

\newcommand*{\rhok}{\rho\iterDown}

\newcommand*{\linFeasLetter}{L}
\newcommand*{\linFeas}{\mathcal{\linFeasLetter}}
\newcommand*{\linFeasApprox}{\mathcal{\linFeasLetter}}
\newcommand*{\linFeasApproxk}{\mathcal{\linFeasLetter}\iterSub}

\newcommand*{\omegak}{\surrogateSign{\omega}\KK}

\newcommand*{\chik}{\chi\KK}
\newcommand*{\chikPre}{\bar\chi\KK}
\newcommand*{\chikNPre}{\bar\chi\KKn}

\newcommand*{\Deltak}{\Delta\itDeltaCmd{\kk}}
\newcommand*{\DeltaMax}{\Delta_{\max}}
\newcommand*{\DeltaMin}{\Delta_{\min}}

\newcommand{\DeltaPre}{\bar{\Delta}}
\newcommand{\DeltaPrek}{\DeltaPre\itDeltaCmd{\kk}}

\newcommand*{\ballLetter}{B}
\newcommand*{\ball}[2]{\ballLetter\args{#1;#2}}
\newcommand*{\trk}{\ballLetter\KK}

\newcommand*{\alphak}{\alpha\iterSub}

\newcommand*{\thetak}{\theta\itthetaCmd\kk}

\newcommand*{\covering}{\mathcal C(\mathcal X)}
\newcommand{\filter}{\mathcal F}
\newcommand*{\filterConst}{\gamma_\theta}

\newcommand*{\CCDelta}{\CC_{\Delta}}
\newcommand*{\CCmu}{\CC_{\mu}}

\newcommand*{\CCdeltan}{\delta_{\mathrm{n}}}
\newcommand*{\CCubn}{\CC_{\mathrm{ubn}}}

\newcommand*{\CCfl}{\mathtt{e}}
\newcommand*{\CCflDiff}{\dot{\CCfl}}
\newcommand*{\LipschitzBoundModel}[1]{L^{#1}}
\newcommand*{\LipschitzBoundModels}{L^{\mfL}}

\newcommand*{\CCflF}{\CCfl_{\fL}}
\newcommand*{\CCflG}{\CCfl_{\gL}}
\newcommand*{\CCflH}{\CCfl_{\hL}}
\newcommand*{\CCflDiffF}{\CCflDiff_{\fL}}
\newcommand*{\CCflDiffG}{\CCflDiff_{\gL}}
\newcommand*{\CCflDiffH}{\CCflDiff_{\hL}}
\newcommand*{\LipschitzBoundModelsF}{\LipschitzBoundModels_{\fL}}
\newcommand*{\LipschitzBoundModelsG}{\LipschitzBoundModels_{\gL}}
\newcommand*{\LipschitzBoundModelsH}{\LipschitzBoundModels_{\hL}}

\newcommand*{\CCnorm}{\mathsf{c}}
\newcommand*{\CCcrit}{\mathtt{w}}

\newcommand*{\armijoConstant}{\mathsf{a}}
\newcommand*{\backtrackConstant}{\mathsf{b}}
\newcommand*{\CCsd}{\mathtt{c}_{\mathrm{sd}}}
\newcommand*{\CCsdT}{\tilde{\mathtt{c}}_{\mathrm{sd}}}

\newcommand*{\CChessF}{\mathtt{H}}

\newcommand*{\chikl}{\chi\itCmd{\kk_{\ell}}}
\newcommand*{\thetakl}{\theta\itthetaCmd{\kk_{\ell}}}
\newcommand*{\vnkl}{\vn\itCmd{\kk_{\ell}}}
\newcommand*{\vskl}{\vs\itCmd{\kk_{\ell}}}
\newcommand*{\vtkl}{\vt\itCmd{\kk_{\ell}}}
\newcommand*{\vxkl}{\vx\itCmd{\kk_{\ell}}}
\newcommand*{\Deltakl}{\Delta\itDeltaCmd{\kk_{\ell}}}
\newcommand*{\vxNkl}{\vxN\itCmd{\kk_{\ell}}}
\newcommand*{\vxSkl}{\vxS\itCmd{\kk_{\ell}}}
\newcommand*{\scalkl}{\scalShort\itCmd{\kk_{\ell}}}
\newcommand*{\mfkl}{\mf\itCmd{\kk_{\ell}}}
\newcommand*{\mgkl}{\mg\itCmd{\kk_{\ell}}}
\newcommand*{\mhkl}{\mh\itCmd{\kk_{\ell}}}
\newcommand*{\linFeasApproxkl}{\mathcal{\linFeasLetter}\iterSubCmd{\kk_{\ell}}}

\newcommand*{\DeltaInit}{\Delta\itDeltaCmd{0}}
\newcommand*{\gammass}{\gamma_{0}}
\newcommand*{\gammas}{\gamma_{1}}
\newcommand*{\gammag}{\gamma_{2}}
\newcommand*{\CCaccept}{\nu_0}
\newcommand*{\CCsuccess}{\nu_1}
\newcommand*{\CCcritTestchi}{\varepsilon_\chi}
\newcommand*{\CCcritTesttheta}{\varepsilon_\theta}

\newcommand*{\CCMu}{\mathtt{M}}
\newcommand*{\CCBeta}{\mathtt{B}}
\newcommand*{\CCalpha}{\alpha}

\newcommand*{\CCcritDenom}{\mathtt{W}}

\newcommand{\vxNext}{\vx^{(\kk+1)}}
%\newcommand{\DeltaNext}{\Delta{(\kk+1)}}
\newcommand{\DeltaPreNext}{\DeltaPre_{(\kk+1)}}
\newcommand{\None}{\mathtt{undef}}

\newcommand*{\acceptIndices}{\mathcal{A}}
\newcommand*{\restoreIndices}{\mathcal{R}}
\newcommand*{\filterIndices}{\mathcal{Z}}

\newcommand*{\CCubj}{\CC_{\mathrm{ubj}}}

\newcommand*{\CCaccN}{\mathtt{e}_{\mathrm{err}}}
\newcommand*{\CCaccS}{\CCaccN}

\newcommand*{\CCthetaub}{\CC_{\mathrm{ub}\theta}}

\newcommand*{\CCepschi}{\epsilon}

\newcommand*{\modd}{\ve{\surrogateSign{d}}}

\newcommand*{\vw}{\ve{w}}

\newcommand*{\front}{\mathcal Y}

\newcommand*{\spreadGamma}{\Gamma}
\newcommand*{\spreadDelta}{\Delta}

\DeclareAcronym{mop}{ short=MOP, long=multi-objective optimization problem}

\DeclareAcronym{moo}{ short=MOO, long=multi-objective optimization}

\DeclareAcronym{rbf}{ short=RBF, long=radial basis function}

\DeclareAcronym{wlog}{ short=w.l.o.g., long=without loss of generality}

\DeclareAcronym{lhs}{ short=LHS, long=left-hand side}

\DeclareAcronym{rhs}{ short=RHS, long=right-hand side}

\DeclareAcronym{mfcq}{ short=MFCQ, long={{M}angasarian-{F}romovitz Constraint Qualification} }

%Better widebar

\makeatletter
\let\save@mathaccent\mathaccent
\newcommand*\if@single[3]{%
  \setbox0\hbox{${\mathaccent"0362{#1}}^H$}%
  \setbox2\hbox{${\mathaccent"0362{\kern0pt#1}}^H$}%
  \ifdim\ht0=\ht2 #3\else #2\fi
  }
%The bar will be moved to the right by a half of \macc@kerna, which is computed by amsmath:
\newcommand*\rel@kern[1]{\kern#1\dimexpr\macc@kerna}
%If there's a superscript following the bar, then no negative kern may follow the bar;
%an additional {} makes sure that the superscript is high enough in this case:
\newcommand*\widebar[1]{\@ifnextchar^{{\wide@bar{#1}{0}}}{\wide@bar{#1}{1}}}
%Use a separate algorithm for single symbols:
\newcommand*\wide@bar[2]{\if@single{#1}{\wide@bar@{#1}{#2}{1}}{\wide@bar@{#1}{#2}{2}}}
\newcommand*\wide@bar@[3]{%
  \begingroup
  \def\mathaccent##1##2{%
%Enable nesting of accents:
    \let\mathaccent\save@mathaccent
%If there's more than a single symbol, use the first character instead (see below):
    \if#32 \let\macc@nucleus\first@char \fi
%Determine the italic correction:
    \setbox\z@\hbox{$\macc@style{\macc@nucleus}_{}$}%
    \setbox\tw@\hbox{$\macc@style{\macc@nucleus}{}_{}$}%
    \dimen@\wd\tw@
    \advance\dimen@-\wd\z@
%Now \dimen@ is the italic correction of the symbol.
    \divide\dimen@ 3
    \@tempdima\wd\tw@
    \advance\@tempdima-\scriptspace
%Now \@tempdima is the width of the symbol.
    \divide\@tempdima 10
    \advance\dimen@-\@tempdima
%Now \dimen@ = (italic correction / 3) - (Breite / 10)
    \ifdim\dimen@>\z@ \dimen@0pt\fi
%The bar will be shortened in the case \dimen@<0 !
    \rel@kern{0.6}\kern-\dimen@
    \if#31
      \overline{\rel@kern{-0.6}\kern\dimen@\macc@nucleus\rel@kern{0.4}\kern\dimen@}%
      \advance\dimen@0.4\dimexpr\macc@kerna
%Place the combined final kern (-\dimen@) if it is >0 or if a superscript follows:
      \let\final@kern#2%
      \ifdim\dimen@<\z@ \let\final@kern1\fi
      \if\final@kern1 \kern-\dimen@\fi
    \else
      \overline{\rel@kern{-0.6}\kern\dimen@#1}%
    \fi
  }%
  \macc@depth\@ne
  \let\math@bgroup\@empty \let\math@egroup\macc@set@skewchar
  \mathsurround\z@ \frozen@everymath{\mathgroup\macc@group\relax}%
  \macc@set@skewchar\relax
  \let\mathaccentV\macc@nested@a
%The following initialises \macc@kerna and calls \mathaccent:
  \if#31
    \macc@nested@a\relax111{#1}%
  \else
%If the argument consists of more than one symbol, and if the first token is
%a letter, use that letter for the computations:
    \def\gobble@till@marker##1\endmarker{}%
    \futurelet\first@char\gobble@till@marker#1\endmarker
    \ifcat\noexpand\first@char A\else
      \def\first@char{}%
    \fi
    \macc@nested@a\relax111{\first@char}%
  \fi
  \endgroup
}
\makeatother

\let\oldbar\bar
\renewcommand*{\bar}[1]{{\mathchoice{\widebar{#1}}{\widebar{#1}}{\widebar{#1}}{\oldbar{#1}}}}

\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}}
\newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}}
\newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}}
\newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}}
\let\oldtilde\tilde
\renewcommand*{\tilde}[1]{\mathchoice{\widetilde{#1}}{\widetilde{#1}}{\oldtilde{#1}}{\oldtilde{#1}}}
\let\oldhat\hat
\renewcommand*{\hat}[1]{\mathchoice{\widehat{#1}}{\widehat{#1}}{\oldhat{#1}}{\oldhat{#1}}}
\let\oldforall\forall
\renewcommand*{\forall}{\mathrel{\oldforall}}


\title[Filter Trust-Region Algorithm for MOPs with Constraints]{Multi-Objective Trust-Region Filter Method for Nonlinear Constraints using Inexact Gradients}

\author{\firstname{Manuel} \lastname{Berkemeier}}
\address{ Department of Computer Science\\
TU Dortmund University\\
Germany}

\thanks{Research funded by the DFG Priority Programme 2353 “Daring more intelligence” (project ID 501834605). We warmly thank the two anonymous reviewers for their insightful input, which has helped to vastly improve this article.}
\CDRGrant[DFG]{501834605}

\email{manuel.berkemeier@tu-dortmund.de}
\author{\firstname{Sebastian} \lastname{Peitz}}
\address{ Department of Computer Science\\
TU Dortmund University\\
Germany}
\email{sebastian.peitz@tu-dortmund.de}


\keywords{\kwd{Multi-Objective Optimization}
\kwd{Multiobjective Optimization}
\kwd{Nonlinear Optimization}
\kwd{Derivative-Free Optimization}
\kwd{Trust-Region Method}
\kwd{Surrogate Models}
\kwd{Filter Method}}


\begin{abstract}
In this article, we build on previous work to present an optimization algorithm for non-linearly constrained multi-objective optimization problems. The algorithm combines a surrogate-assisted derivative-free trust-region approach with the filter method known from single-objective optimization. Instead of the true objective and constraint functions, so-called \emph{fully linear} models are employed, and we show how to deal with the gradient inexactness in the composite step setting, adapted from single-objective optimization as well. Under standard assumptions, we prove convergence of a subset of iterates to a quasi-stationary point and, if constraint qualifications hold, then the limit point is also a KKT-point of the multi-objective problem.
\end{abstract}


\datepublished{2025-12-08}
\begin{document}
\theoremstyle{plain}
\newtheorem{assu}[theo]{Assumption}

\maketitle


\section{Introduction}

\Acfp{mop} arise naturally in many areas of mathematics, engineering, in the natural sciences or in economics. The goal of \acf{moo} is to find acceptable trade-offs between the competing objectives of a
\ac{mop}. Generally, there are multiple solutions constituting the so-called Pareto Set in variable space and the Pareto Front in objective space. In this article, we consider the problem
\begin{equation}
\min_{\vx\in \mathbb{R}^\varDim}
\begin{bmatrix}
\fL_1(\vx)
\\
\vdots
\\
\fL_\fDim(\vx)
\end{bmatrix} =
\min_{\vx \in \mathbb{R}^\varDim} \vf(\vx)
\quad\text{s.t.}\quad
\vh(\vx) = \tran{[\hL_1 (\vx), \ldots, \hL_\hDim(\vx)]} = \vZ,
\;
\vg(\vx) = \tran{[\gL_1 (\vx), \ldots, \gL_\gDim(\vx)]} \le \vZ,
\tag{MOP}\label{eqn:mop_nonlin}
\end{equation}
where all functions are twice continuously differentiable. A global Pareto-optimal point $\vx^*$ is feasible and non-dominated, i.e., there is no other feasible $\vx \ne \vx^*$ with $\vf(\vx) \le \vf(\vx^*)$ and $\fL_\indf(\vx) < \fL_\indf(\vx^*)$ for some $\indf\in\{1,\ldots,\fDim\}$. Throughout this article we will refer to the feasible set as $\feas :=
\set{
\vx \in \mathbb{R}^\varDim :
\vh(\vx) = \vZ \text{ and }
\vg(\vx) \le \vZ }$.

There is a multitude of methods available to approximate a single solution or the entire Pareto Set/Front (or a superset thereof) and the choice of method is heavily dependent on the structure of the problem at hand and the demands of the person seeking a solution. The references~\cite{ehrgottMulticriteriaOptimization2005, eichfelderTwentyYearsContinuous2021a, jahnVectorOptimizationTheory2011, miettinenNonlinearMultiobjectiveOptimization2013} all provide an extensive overview of the topic. Amongst others, there are scalarization approaches~\cite{eichfelderAdaptiveScalarizationMethods2008} and adaptions of single-objective descent methods to the multi-objective case for different problem classes as defined by the properties of their objectives and constraint functions, see~\cite{drummondSteepestDescentMethod2005, fliegeSteepestDescentMethods2000, fliegeMethodConstrainedMultiobjective2016, fukudaSURVEYMULTIOBJECTIVEDESCENT2014, gebkenEfficientDescentMethod2021, gebkenDescentMethodEquality2018, lucambioperezNonlinearConjugateGradient2018, wilppuNewMultipleSubgradient2014} for examples of well-known scalar techniques adjusted for constrained and unconstrained \acp{mop} in the smooth and non-smooth case. For global approximations there are, e.g., evolutionary algorithms (see~\cite{coelloEvolutionaryAlgorithmsSolving2013} for an overview and~\cite{debFastElitistMultiobjective2002} for a prominent example) and structure-exploiting methods~\cite{gebkenHierarchicalStructurePareto2019, hillermeierNonlinearMultiobjectiveOptimization2001, martinParetoTracerPredictor2018}. Of course, there is also research in combining global and local techniques~\cite{peitzGradientBasedMultiobjectiveOptimization2018, schutzeParetoExplorerGlobal2020}. Whilst there are natural applications of \ac{moo} for machine learning tasks~\cite{albuquerqueMultiobjectiveTrainingGenerative2019, senerMultiTaskLearningMultiObjective2019}, machine learning techniques can conversely be used to assist the search for optimal points~\cite{navonLearningParetoFront2021}. In case of expensive objective or constraints, surrogate models can be employed~\cite{chughSurveyHandlingComputationally2019, debSurrogateModelingApproaches2020, peitzSurveyRecentTrends2018}.

In our setting we assume (some) objective and constraint functions to be computationally expensive and without exact derivative information available. This motivates the use of \emph{derivative-free} optimization methods, which also have been adapted to the multi-objective case. For example, there are direct search algorithms~\cite{audetMeshAdaptiveDirect2010, bigeonDMultiMADSMeshAdaptive2021, bigeonHandlingConstraintsMultiobjective2024, custodioDirectMultisearchMultiobjective2011, dedonckerAdaptiveDirectMultisearch2022, silvaInexactRestorationDirect2024}, implicit filtering~\cite{cocchiImplicitFilteringAlgorithm2018} and linesearch~\cite{liuzziDerivativeFreeApproachConstrained2016} algorithms, or trust-region algorithms using derivative-informed approximations~\cite{mohammadiTrustregionApproachComputing2024, quTrustRegionMethods2013, villacortaTrustRegionMethodUnconstrained2014} or derivative-free surrogate models~\cite{mohammadiTrustregionApproachComputing2025, ryuDerivativeFreeTrustRegionMethod2014, thomannTrustRegionAlgorithmHeterogeneous2019}. Based on such trust-region algorithms, we have in a previous article~\cite{berkemeierDerivativeFreeMultiobjectiveTrust2021} presented a trust-region algorithm for problems with a feasible set that is convex and compact (or $\feas = \mathbb{R}^\varDim$). The algorithm uses fully linear models (e.g., Lagrange polynomials or \acp{rbf} as in~\cite{wildORBITOptimizationRadial2008}) to approximate the objective functions and passes the exact constraints to an inner solver. We build upon this work to accommodate general non-linear constraints by also modelling them and solving inexact sub-problems.

To this end, we transfer the techniques from~\cite{fletcherGlobalConvergenceTrustRegion2002} to the multi-objective case with inexact derivatives. Inexact derivative have already been handled in single-objective optimization in a similar manner. For example, the authors of~\cite{ferreiraGlobalConvergenceDerivativefree2017} provide strong convergence results for inexact objectives and objective derivatives with similar model accuracy requirements to our case. We will also discuss the single-objective algorithm presented in~\cite{easonTrustRegionFilter2016}, as it is also based on fully linear models and employs a special kind of criticality check. Our work, however, is more along the lines of~\cite{waltherInexactTrustregionSQPfilter2016}, where the filter trust-region algorithm and the composite step framework are likewise modified to handle inexact derivatives of both the objective and constraint functions. But in~\cite{waltherInexactTrustregionSQPfilter2016} the derivatives (approximated via automatic differentiation) can eventually become exact, in contrast to our setting, which is why (without constraint qualifications) we can only prove convergence to a \emph{quasi-stationary} point, similar to~\cite{echebestInexactRestorationDerivativefree2017}. The algorithms in~\cite{echebestInexactRestorationDerivativefree2017} and~\cite[Chapter~4]{garreisOptimalControlUncertainty2019} also use fully linear surrogate models satisfying the same error bounds as assumed in this article within single-objective filter algorithms. In~\cite{echebestInexactRestorationDerivativefree2017}, convergence to KKT-points is proven under Linear Independence Constraint Qualifications. Similarly, we will see that our convergence result can be proven under the \aclp*{mfcq}. Such constraint qualifications are not needed in~\cite{garreisOptimalControlUncertainty2019}. But to show convergence to critical points, in~\cite{garreisOptimalControlUncertainty2019} an additional bound for the model accuracy is required, which is viable if the projection onto the feasible set can be done with arbitrary precision. In contrast, our approach uses projections on approximations of the linearized feasible set.

In next section, we will state the relevant optimality conditions for the multi-objective case. Afterwards, the main building blocks of our algorithm are described in~\cref{sec:trm_ideas}. The algorithm itself is given in~\cref{section:algorithm} and convergence is shown in~\cref{sec:quasi_convergence,section:kkt_convergence}. Finally, numerical examples and comparisons are presented in~\cref{sec:examples} with a brief discussion in~\cref{sec:conclusion}.

\subsection{Optimality Conditions}

In this subsection we introduce necessary conditions for a point to be locally Pareto-optimal according to the following definition.

\begin{defi}\label{def:locally_optimal}
The point $\vx^*\!\in\Mk\feas$ is locally Pareto-optimal, if there is a neighborhood $\mathcal U_\varepsilon = \{\vx \Mk\in\Mk \feas: \norm{\vx \mk-\mk \vx^*} \Mk\leq\Mk \varepsilon\}$ such that $\vx^*$ is non-dominated in $\mathcal U_\varepsilon$, that is, there is no $\vx \in \mathcal U_\varepsilon$ such that for all $\indf \in \{1, \ldots, \fDim\}$ it holds that $\fL_\indf(\vx) \leq \fL_\indf(\vx^*)$ and $f_{\indmax{\indf}}(\vx) < f_{\indmax{\indf}}(\vx^*)$ for at least one index $\indmax{\indf} \in \{1, \ldots, \fDim\}$.
\end{defi}

To state the optimality criteria, we require the following assumption to hold, so that all functions are sufficiently smooth:

\begin{assu}\label{ass:lipschitz_gradients_exact}	
The objective functions $\fL_\indf,\, \indf=1,\ldots,\fDim$, and the constraint functions $\gL_\indg,\, \indg=1, \ldots, \gDim$, and $\hL_\indh,\,\indh=1,\ldots,\hDim$, are twice continuously differentiable in an open domain containing $\feas$ and have Lipschitz continuous gradients on $\feas$.
\end{assu}

There then is a formulation of Fermat's Theorem for optimization with multiple objectives:

\begin{theo}[{\cite[Theorem~5.2]{linMaximalVectorsMultiobjective1976}}]\label{thm:fermats_theorem}
Suppose \cref{ass:lipschitz_gradients_exact} holds. If $\vx^* \in \feas$ is locally Pareto-optimal for~\eqref{eqn:mop_nonlin}, then there is no $\vd \in \tangcone{\feas}{\vx^*}$ such that for all $\indf\in \{1,\ldots, \fDim\}$ it holds that
\(
\dotProd{\vd}{\df_\indf(\vx^*)} < 0,
\) where $\tangcone{\feas}{\vx^*}$ is the tangent cone of $\feas$ at $\vx^*$, i.e., the set of vectors $\ve{v} \in \mathbb{R}^\varDim$ for which there is a sequence $\vvk\in\feas$ and scalars $\alphak > 0$ with
\[
\lim_{\kk\to\infty} \vvk = \vx^*,
\lim_{\kk\to\infty} \alphak = 0,
\lim_{\kk\to\infty} \frac{\vvk - \vx^*}{\alphak} = \ve{v}.
\]
\end{theo}

We call a point $\vx^* \in \feas$ satisfying the criterion in \cref{thm:fermats_theorem} \emph{Pareto-critical} or stationary. \Cref{thm:fermats_theorem} can be used to motivate the following problem to compute a descent direction and check for criticality:
\begin{equation}
\min_{\vd \in \tangcone{\feas}{\vx}, \normTwo{\vd} \leq 1}
\max_{\indf = 1, \ldots, \fDim}
\dotProd{\vd}{\df_\indf (\vx)}\label{eqn:DemandIdenticalOval}
\end{equation}
For a critical point $\vx\in\feas$ the optimal value is zero, else the minimizer is a multi-descent direction (cf.~\cite[Theorem~1.8.]{attouchDynamicGradientApproach2015}). The choice of norm in problem~\eqref{eqn:DemandIdenticalOval} is not really important. We could instead use a linear norm and assume that constraint qualifications hold at $\vx \in \feas$, ensuring that the tangent cone equals the set of linearized directions
\[
L(\vx) =
\set{
\vd \in \mathbb{R}^\varDim :
\dotProd{\vd}{\dh_\indh(\vx)} = 0, \, \indh = 1, \ldots, \hDim,
\;
\dotProd{\vd}{\dg_\indg(\vx}) \le 0, \, \indg \in A(\vx) },
\]
with $A(\vx) := \{\indg \in \{1,\ldots,\gDim\}: g_\indg(\vx) = 0 \}$, to obtain a linear problem. These constraint qualifications are also used to derive the KKT conditions for~\eqref{eqn:mop_nonlin} from~\cref{thm:fermats_theorem}~\cite[Theorem~2.3.]{marusciaciFritzJohnType1982}. See~\cite{fliegeSteepestDescentMethods2000, mukaiAlgorithmsMulticriterionOptimization1980} for prominent examples of descent algorithms using the set of linearized directions. In our algorithm, however, we do not use the sets $L(\vxk)$. Instead, at $\vxk$ (which does not have to be feasible) we use the linearized feasible set,
\[
\linFeas(\vxk) :=
\set{
\vxk + \vd \in \mathbb{R}^\varDim:
\vh(\vxk) + \dotMatProd{\DH(\vxk)}{\vd} = \vZ, \;
\vg(\vxk) + \dotMatProd{\DG(\vxk)}{\vd} \le \vZ },
\]
where $\DH(\vxk)$ and $\DG(\vxk)$ denote the \emph{full} Jacobian matrices of the constraints $\vh$ and $\vg$, respectively. The set $\linFeas(\vxk) - \vxk$ is not necessarily a cone, but for feasible $\vxk$, it is a subset of $L(\vxk)$, and intuitively it should not matter which of the sets is used near critical points. Indeed, the set $\linFeas(\vxk)$ appears in the sequential optimality condition in~\cite{martinezPracticalOptimalityCondition2003}. Similar to~\cite{echebestInexactRestorationDerivativefree2017}, we are hindered from using this sequential condition directly by the fact that $\vg$ and $\vh$ are actually approximated by local surrogate models (see next section). For now, assume that $\linFeasApproxk$ is a surrogate-based \emph{approximation} of $\linFeas(\vxk)$. We will show that our algorithm produces a sequence that approximates a feasible point, which is also \emph{quasi-stationary}: the values of~\eqref{eqn:DemandIdenticalOval} tend to zero if we replace $\tangcone{\feas}{\vxk}$ with $\linFeasApproxk$. The convergence result is strengthened by requiring the \acp{mfcq} to hold at the limit point $\bar{\vx}$. Then the solution set sequence of the problems resulting from substituting $\linFeasApproxk$ into~\eqref{eqn:DemandIdenticalOval} is continuous at $\bar{\vx}$, and a vanishing approximation error for $\linFeasApproxk$ ensures that in the limit the solution set is the same as it is for the exact linearized set $\linFeas(\bar{\vx})$. Furthermore, the \acp{mfcq} make the limit point a KKT point (similar to~\cite{gonzagaGloballyConvergentFilter2004}).

\begin{theo}\label{thm:subproblems_implies_kkt}
Suppose that $\vx$ is feasible and that suitable constraint qualifications hold, e.g., the \acp{mfcq}. If the linear optimization problem
\begin{equation}
-
\min_{
\substack{
\vd \in \linFeas(\vx) - \vx,\\
\normInf{\vd} \leq 1 } }
\max_{\indf = 1, \ldots, \fDim}
\dotProd{\vd}{\df_\indf (\vx)}\label{eqn:linear_opt_problem}
\end{equation}
has zero as its optimal value, then $\vx$ is also a KKT-point of~\eqref{eqn:mop_nonlin}.
\end{theo}

\begin{proof}
We refrain from listing the full set of KKT conditions and instead refer to
\begin{itemize}
\item their first appearance for \acp{moo} with inequality constraints,~\cite{kuhnNonlinearProgramming2014},
\item an extension to equality constraints,~\cite{marusciaciFritzJohnType1982},
\item and~\cite{hillermeierNonlinearMultiobjectiveOptimization2001, martinParetoTracerPredictor2018} for more contemporary notation.
\end{itemize}

Dropping the argument $\vx$ for notational convenience, e.g., letting $\vg = \vg(\vx) \leq 0$ etc., and denoting by $\DF,\DH, \DG$ the objective and constraint Jacobians, the linear problem~\eqref{eqn:linear_opt_problem} is equivalent to
\begin{equation}
\max_{\vd, \beta^-}
\begin{bmatrix}
\tran{\vZ_\varDim} & 1
\end{bmatrix}
\begin{bmatrix}
\vd \\ \beta^-
\end{bmatrix}
\quad\text{s.t.}\quad
\vd \in \mathbb{R}^\varDim,
\;
\beta^- \in \mathbb R,
\;
\begin{bmatrix}
-\eye{\varDim} & \vZ_\varDim \\
\eye{\varDim} & \vZ_\varDim \\
\DF & \ve{1}_\fDim \\
\DH & \vZ_\hDim \\
\DG & \vZ_\gDim
\end{bmatrix}
\begin{bmatrix}
\vd \\ \beta^-
\end{bmatrix}
\begin{matrix}
\le \\ \le \\ \le \\ = \\ \le
\end{matrix}
\begin{bmatrix}
\ve{1}_\varDim\\
\ve{1}_\varDim\\
\vZ_\fDim \\
\vZ_\hDim\\
-\vg \\
\end{bmatrix}.\label{eqn:legitimizersSannop}
\tag{P}
\end{equation}
Consider also the dual problem:
\begin{equation}
\begin{aligned}
&\min_{\ve{y}^1, \ldots, \ve{y}^5 }
\begin{bmatrix}
\tran{\ve{1}_{\varDim}} & \tran{\ve{1}}_{\varDim} & \tran{\vZ}_{\fDim} & \tran{\vZ}_{\hDim} & - \tran{\vg}
\end{bmatrix}
\begin{bmatrix}
\ve{y}^1 \\ \vdots \\ \ve{y}^5
\end{bmatrix}
\quad
\text{s.t.}
\quad
\ve{y}^1 \ge \vZ_{\varDim},
\ve{y}^2 \ge \vZ_{\varDim},
\ve{y}^3 \ge \vZ_{\fDim},
\ve{y}^4 \in \mathbb{R}^\hDim,
\ve{y}^5 \ge \vZ_{\gDim}
\quad
\text{and}
\\
&
\begin{bmatrix}
-\eye{\varDim} & \eye{\varDim} & \tran{\DF} & \tran{\DH} & \tran{\DG}\\
\tran{\vZ}_{\varDim} & \tran{\vZ}_{\varDim} & \tran{\ve{1}}_{\fDim} & \tran{\vZ}_{\hDim} & \tran{\vZ}_{\gDim}
\end{bmatrix}
\begin{bmatrix}
\ve{y}^1 \\ \vdots \\ \ve{y}^5
\end{bmatrix} =
\begin{bmatrix}
\vZ_{\varDim} \\ 1
\end{bmatrix}.
\end{aligned}
\tag{D}\label{eqn:buntingGeologically}
\end{equation}
If $0$ is the optimal value of~\eqref{eqn:legitimizersSannop}, then $\beta^-=0$. By strong duality, the dual problem is feasible with optimal value 0, implying $\ve{y}^1 = \vZ$, $\ve{y}^2 = \vZ$ and $-\noDotProd{\ve{g}}{\ve{y}^5}= 0$. The KKT equations immediately follow from the remaining constraints and from the complementary slackness property of dual solution pairs.
\end{proof}

\section{Trust-Region Concepts and Surrogates}\label{sec:trm_ideas}
As mentioned in the previous section, we assume (at least some) of the objectives or constraints to be computationally expensive. To approximate a Pareto-critical point whilst avoiding expensive function evaluations, a trust-region approach is used: The true functions $\vf,\vh$ and $\vg$ are modeled by $\mfk,\mhk$ and $\mgk$ respectively. These models are constructed to be sufficiently accurate within iteration-dependent trust-regions
\begin{equation}
\trk = B(\vxk; \Deltak) :=
\set{
\vx \in \mathbb{R}^\varDim : \normTrIt{ \vx - \vxk } \le \Deltak}.\label{eqn:bonhamPumpkin}
\end{equation}

Using the surrogate models, a step $\vsk$ and a step-size $\stepsizek$ are determined, in such a way as to reduce the constraint violation or achieve an objective value reduction at the trial point $\vxk + \stepsizek \vsk$, which is tested as a candidate for the next iterate. The step $\stepsizek\vsk$ is computed so that the trial point is contained in both the trust-region and the approximate linearized feasible set at $\vxk$:
\begin{equation}
\linFeasApproxk =
\linFeasApproxk(\vxk) =
\set{
\vxk + \vs \in \mathbb{R}^\varDim:
\mhk(\vxk) + \DmH(\vxk) \cdot \vs = \vZ, \;
\mgk(\vxk) + \DmG(\vxk)\cdot \vs \le \vZ},\label{eqn:orgiastThripses}
\end{equation}
where $\DmH$ and $\DmG$ are now the full Jacobians of the constraint function models. This introduces uncertainty and iterates might no longer be feasible. Hence, we treat the constraints as \emph{relaxable}, and we might have to evaluate $\vf$ (and $\vg$ and $\vh$) outside of $\feas$, which motivates the next assumption. It bounds the trust-regions and ensures that all functions are available for all possible iterates, trust-region norms and radii. Note, that the algorithm actually can use two trust-region sizes, the preliminary size $\DeltaPrek$ at the beginning of an iteration and $\Deltak \leq \DeltaPrek$ after the \emph{Criticality Routine}, which is described in \cref{section:crit_routine}.

\begin{assu}\label{ass:covering}
There is a constant $0<\DeltaMax<\infty$ such that for every $\kk \in \mathbb{N}_0$ the trust region sizes conform to $0 < \Deltak \le \DeltaPrek < \DeltaMax$. All functions (the true functions and their models) are defined in these regions, i.e., on the set $\covering = \bigcup_{\kk\in \mathbb{N}_0}
\ball{\vxk}{\DeltaPrek} $.
\end{assu}

Note, that the trust-regions in~\eqref{eqn:bonhamPumpkin} may be defined using iteration-dependent norms $\normTrIt{\bullet}$. Later, we will also introduce norm-constrained subproblems, where the norm $\normIt{\bullet}$ can vary, too. This allows for great flexibility in the actual implementation. We only require that the norms be uniformly equivalent to the Euclidean norm:

\begin{assu}\label{ass:equivalent_norms}
There is a constant $\CCnorm \ge 1$ such that for all $\vx\in \covering$ and all $\kk\in \mathbb{N}_0$ and $\norm{\bullet}_* = \normTrIt{\bullet}$ or $\norm{\bullet}_* = \normIt{\bullet}$ it holds that
\begin{equation}
\frac{1}{\CCnorm}
\norm{\vx}_*
\le
\normTwo{\vx}
\le
\CCnorm \norm{\vx}_*.\label{eqn:equivalent_norms}
\end{equation}
\end{assu}

\begin{rema}\label{remark:equivalent_norms}
Any two norms that are uniformly equivalent to $\normTwo{\bullet}$ with constant $\CCnorm$ are pairwise equivalent with constants $\nicefrac{1}{\CCnorm^2}$ and $\CCnorm$.
\end{rema}

\subsection{Fully Linear Models}

In this subsection, we first want to explain what it means for the model functions to be sufficiently accurate. Although the true derivative information is not used in any computation, we will need it for the convergence analysis and hence require the following generalization of~\cref{ass:lipschitz_gradients_exact}:

\begin{assu}\label{ass:lipschitz_gradients_true}
The objective functions $\fL_\indf,\, \indf=1,\ldots,\fDim$, the constraint functions $\gL_\indg,\, \indg=1, \ldots, \gDim$, and $\hL_\indh,\,\indh=1,\ldots,\hDim$, are twice continuously differentiable in an open domain containing $\covering$ and have Lipschitz continuous gradients on $\covering$.
\end{assu}

With \cref{ass:lipschitz_gradients_true} we can use~\cite[Definition~6.1.]{connIntroductionDerivativefreeOptimization2009} to have our models satisfy error bounds w.r.t. to their construction radius:

\begin{defi}[Fully Linear Models]\label{def:fully_linear_models}
Let $\DeltaMax > 0$ be a given constant and let $\fL\colon \mathbb{R}^n\to \mathbb{R}$ be a scalar-valued function satisfying \cref{ass:lipschitz_gradients_true}. A set of model functions $\mathcal M = \{ \mfL \colon \mathbb{R}^n \to \mathbb{R} \}$ is called a fully linear class of models for $f$ if the following hold:
\begin{enumerate}
\item\label{def:fully_linear_models1}
There are positive constants $\CCfl,\CCflDiff$ and $\LipschitzBoundModels$ such that for any $\Delta\in (0,\DeltaMax]$ and for any $\vx \in \covering$ there is a model function $\mfL \in \mathcal M$, with Lipschitz continuous gradient and corresponding Lipschitz constant bounded by $\LipschitzBoundModels$, such that the error between the gradients and the error between the values satisfy
\begin{equation*}
\normTwo{\df(\ve {\xi}) - \dmf(\ve {\xi})}
\le \CCflDiff \Delta
\quad\text{and}\quad
\abs{
\fL(\ve {\xi}) - \mfL(\ve {\xi})}
\le \CCfl \Delta^2, \qquad\forall\ve {\xi}
\in \ball{\vx}{\Delta} \cap \covering.
\end{equation*}
\item\label{def:fully_linear_models2}
For this class $\mathcal M$ there exists a ``model-improvement'' algorithm that, in a finite, uniformly bounded number of steps can either establish that a given model $\mfL\in \mathcal M$ is fully linear on $\ball{\vx}{\Delta}$ or find a model $\mfL \in \mathcal M$ that is fully linear on $\ball{\vx}{\Delta}$.
\end{enumerate}
\end{defi}

\begin{defi}\label{def:fully_linear_vector_models}
Let $\DeltaMax>0$ be a given constant and let $\vf = \tran{[\fL_1, \ldots, \fL_\fDim]}$ be a vector of $\fDim\geq 1$ functions satisfying the requirements of
\cref{def:fully_linear_models} with classes $\mathcal M_\indf$ and constants $(\CCfl_\indf, \CCflDiff_\indf, \LipschitzBoundModels_\indf), \indf=1,\ldots,\fDim$. Then
\[
\mathcal M =
\{ \mf = \tran{[\mfL_1, \ldots, \mfL_\fDim]} :
\mfL_1 \in \mathcal M_1, \
\ldots, \
\mfL_\fDim \in \mathcal M_\fDim \}
\]
is a class of fully linear, vector-valued model functions with constants $\max_\indf \CCfl_\indf > 0$, $\max_\indf \CCflDiff_\indf>0$ and $\max_\indf\LipschitzBoundModels_\indf>0$. A vector of functions $\mf \in \mathcal M$ is deemed fully linear, if all components $\mfL_\indf$ are fully linear. The improvement algorithms for $\mathcal M_\indf$ are applied component-wise.
\end{defi}

\begin{rema}\label{rem:radius_dependent_surrogates}
In trust-region algorithms with derivatives, oftentimes Taylor polynomial approximations are used, and under standard assumptions they satisfy the error bounds. In fact, the error bounds are motivated by these traditional models. In contrast, the notion of fully linear models allows for more flexible surrogates, the construction of which usually depends not only on the iterate but also the current trust region radius.
\end{rema}

A function that is fully linear for $\Deltak > 0$ is automatically fully linear for any smaller radius $0<\Delta<\Deltak$. When the trust-region radius is bounded above, then the constants $\CCfl>0$ and $\CCflDiff > 0$ can be chosen large enough such that a fully linear model stays fully linear in enlarged trust-regions:

\begin{lemm}[{\cite[Lemma~10.25]{connIntroductionDerivativefreeOptimization2009}}]\label{thm:fully_linear_enlarged_trust_region}
For $\vx \in \covering$ and $\Delta \in (0,\DeltaMax]$ consider a function $\fL$ and a fully linear model $\mfL$ of $\fL$ with constants $\CCfl, \CCflDiff, \LipschitzBoundModels > 0$. Let $L^{\fL} > 0$ be a Lipschitz constant of $\df$. Assume \ac{wlog} that
\begin{equation}
\LipschitzBoundModels + L^{\fL} \le \CCfl
\text{ and }
\nicefrac{\CCflDiff}{2} \le \CCfl.\label{eqn:large_fully_linear_constants}
\end{equation}
Then $\mfL$ is fully linear on $\ball{\vx}{\tilde{\Delta}}$ for any $\tilde{\Delta} \in [\Delta, \DeltaMax]$ with the same constants.
\end{lemm}

\begin{assu}\label{ass:fully_linear_models}
For any $\kk \in \mathbb{N}_0$ the models $\mfk$ are fully linear on $\ball{\vxk}{\Deltak}$ as in \cref{def:fully_linear_vector_models} with constants $\CCflF, \CCflDiffF$ and $\LipschitzBoundModel{\vf}$ that are chosen large enough such that~\eqref{eqn:large_fully_linear_constants} is fulfilled globally. The same holds for the models $\mgk$ of $\vg$ with constants $\CCflG, \CCflDiffG, \LipschitzBoundModel{\vg}$ and the models $\mhk$ of $\vh$ with constants $\CCflH, \CCflDiffH$ and $\LipschitzBoundModel{\vh}$. We also assume that all models are \emph{interpolating} at $\vxk$.
\end{assu}

\subsection{Composite Step Approach and Sub-Problems} Just like in many previous articles on multi-objective trust-region algorithms~\cite{berkemeierDerivativeFreeMultiobjectiveTrust2021, quTrustRegionMethods2013, thomannTrustRegionAlgorithmHeterogeneous2019, villacortaTrustRegionMethodUnconstrained2014}, we use the maximum-scalarization
\[
\scalShort(\vx) = \scal{\vf}(\vx) := \max_{\indf=1,\ldots,\fDim} \fL_\indf(\vx)
\quad\text{and}\quad
\scalk(\vx) = \scal{\mfk}(\vx) := \max_{\indf=1,\ldots,\fDim} \mfL_\indf(\vx)
\]
to determine objective reduction and prove convergence. The idea then is to find a step $\vsk$ that approximately solves
\begin{equation*}
\min_{\vs \in \mathbb{R}^\varDim}
\scal{\mfk}(\vxk + \vsk)
\quad\text{s.t.}\quad
\vs \in \args{ \linFeasApproxk - \vxk},
\,
\normTrIt{\vs} \le \Deltak,
\end{equation*}
with the inexact linearized feasible set defined in~\eqref{eqn:orgiastThripses}. Without constraints, inexact line-search can be used. In our case, we use composite-step approach and $\vsk$ is split into a normal component $\vnk$ towards feasibility and a descent direction $\vdk$ (see~\cite{fletcherGlobalConvergenceTrustRegion2002,waltherInexactTrustregionSQPfilter2016} for details). Then, the normal step can be computed with
\begin{equation}
\min \normTwo{ \vn}^2
\quad\text{s.t.}\quad
\vn \in \args{ \linFeasApproxk - \vxk}.\label{eqn:itr_normal}
\tag{ITRN\ensuremath{\KK}}
\end{equation}
If a normal step $\vnk$ has been found (and if $\normTrIt{\vn} \le \DeltaPrek$) the descent direction $\vdk$ can be taken as the minimizer of
\begin{equation}\label{eqn:itr_tangential}
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normIt{\bullet}} =-
\min_{
\beta \in \mathbb{R},
\vd \in \mathbb{R}^\varDim}
\beta
\quad\text{s.t.}\quad
\DmF(\vxNk) \cdot \vd \le \beta, \
\vd \in \args{\linFeasApproxk - \vxNk}, \
\normIt{\vd} \le 1,
\tag{ITRT\ensuremath{\KK}}
\end{equation}
where $\linFeasApproxk$ is the approximate linearized feasible set at $\vxk$ according to~\eqref{eqn:orgiastThripses} and $\vxNk$ is the point $\vxNk = \vxk + \vnk$.
We actually only compute the descent direction if there is enough movement possible after performing the normal step, i.e., the normal step must not be too large. We call such a step \emph{compatible}, and it is defined with respect to the preliminary radius $\DeltaPrek$ :
\begin{defi}\label{def:compatible_normal_step}
Let $\CCDelta \in (0,1], \CCmu > 0$ and $\mu\in (0,1)$ be constants. The minimizer $\vnk$ of~\eqref{eqn:itr_normal} is called compatible (w.r.t. $\DeltaPrek$) if
\begin{equation}
\normTrIt{\vnk}
\le
\CCDelta \DeltaPrek
\min \set{ 1, \CCmu \DeltaPrek^\mu}.\label{eqn:compatible_normal_step}
\end{equation}
Similarly, we call~\eqref{eqn:itr_normal} compatible if it has a non-empty solution-set and a minimizer that is compatible.
\end{defi}

Furthermore, we assume a normal step to exist if the true constraint violation is not too large, as measured by the infeasibility function
\[
\theta\args{\vx} :=
\max\set{ 0,
\max_{\indh=1, \ldots, \hDim}
\abs{ \hL_\indh(\vx)},
\max_{\indg=1,\ldots,\gDim}
\gL_\indg(\vx)}.
\]

\begin{assu}[Existence and Boundedness of Normal Step]\label{ass:normal_step_exists_bounded}
If $\thetak = \theta\args{ \vxk} \le \CCdeltan$, for a constant $\CCdeltan > 0$, then $\vnk$ exists and there is a constant $\CCubn > 0$ such that
\begin{equation}
\normTrIt{\vnk} \le \CCubn \thetak.\label{eqn:normal_step_exists_bounded}
\end{equation}
\end{assu}

\Cref{ass:normal_step_exists_bounded} is a standard assumption and the reasoning behind it can be found in~\cite{fletcherGlobalConvergenceTrustRegion2002}. Finally, the step-size $\stepsizek$ is determined in such a way that $\normTrIt{\vxk + \vnk + \stepsizek \vdk} \le \Deltak$ and a \emph{sufficient decrease condition} (for the objective surrogates) is satisfied, which is described and justified in \cref{section:sufficient_decrease}. The term $\stepsizek \vdk$ is also called \emph{tangential step} and altogether the step $\vsk$ results in the \emph{trial point} $\vxSk = \vxk + \vsk = \vxNk + \stepsizek \vdk$.

\subsection{Filter Mechanism}\label{sec:filter_mechanism}
As mentioned already, iterates can become infeasible. We employ a so-called \emph{filter} to drive them back towards the feasible set. Thus, the trial point $\vxSk = \vxk + \vsk$ is tested not only for actual decrease of the original functions but also against previous iterates stored in the filter $\filter$. $\filter$ is a set of tuples $(\theta_j, \Phi_j), j = 1, \ldots, \abs{\filter}$, describing a forbidden area in image space. In fact, the tuples in $\filter$ (w.r.t. $\vf$) are currently non-dominated for the bi-objective optimization problem of minimizing both $\theta(\vx)$ and $\scalf{\vx}$. The trial point $\vxSk$ is only acceptable for $\filter$ if its value tuple is also (sufficiently) non-dominated. An acceptable trial point is further tested, and if it sufficiently reduces the \emph{true} objectives it is kept as the next iterate. The current iterate might be added to the filter if the predicted objective decrease is small compared to the constraint violation, or if the latter is too large. More formally, we use the following definition:

\begin{defi}[Multi-Objective Filter]\label{def:filter}
A filter $\filter$ with constant $\filterConst \in (0,1)$ for a function $\vf\colon \mathbb{R}^\varDim \to \mathbb{R}^\fDim$ is a discrete set of tuples $\{(\theta_j, \Phi_j)\}_j\subset \mathbb{R}^2$ such that for all $j$ the tuple $(\theta_j, \Phi_j) \in \filter$ is acceptable for $\filter\setminus\set{(\theta_j,\Phi_j)}$. A tuple $(\theta_j, \Phi_j)$ is \emph{acceptable} for $\filter$ iff
\begin{equation*}
\theta_j\le (1-\filterConst) \theta_i
\quad\text{or}\quad
\Phi_j\le \Phi_i - \filterConst \theta_i
\quad \forall (\theta_i, \Phi_i)\in \filter.
\end{equation*}
If a tuple $(\theta_j, \Phi_j)$ is added to $\filter$, then all redundant entries are removed, i.e., all tuples $(\theta_i, \Phi_i)$ with
\begin{equation*}
\theta_i \geq \theta_j
\quad\text{and}\quad
\Phi_i - \filterConst\theta_i
\geq
\Phi_j - \filterConst\theta_j.
\end{equation*}
A \emph{point} $\vx$ is acceptable for the filter $\filter$ iff the value tuple $(\theta\args{\vx}, \scalf{\vx})$ is acceptable for $\filter$, and a point $\vx$ is added to $\filter$ by adding its value tuple.
\end{defi}

As can be seen from~\cref{def:filter}, a filter strengthens non-dominance testing by employing a positive offset $\filterConst \theta_j$. Like with the extended definition of (weak) dominance in~\cite{fliegeMethodConstrainedMultiobjective2016}, feasible and infeasible points can be compared (in contrast to the alternative definition in~\cite{bigeonHandlingConstraintsMultiobjective2024}). Note, that we use a $2$-dimensional filter in accordance with the maximum-scalarization. This makes the convergence analysis easier. It should, however, be possible to use a stricter, $(\fDim + 1)$-dimensional filter (similar to~\cite{gouldMultidimensionalFilterAlgorithm2004}) by using $\vf(\bullet)$ instead of $\scalf{\bullet}$. In~\cref{sec:outerAlgo}, there is a simple draft of such a procedure. Likely, this requires further modifications to the algorithm, e.g., when testing if a step reduces the objective values to a satisfactory degree. For now, we stick to the procedure given in the next section.


\subsection{Sufficient Decrease}\label{section:sufficient_decrease}

In this subsection we want to explain what is meant by ``sufficient decrease''. In short, we have to relate the model predicted objective reduction to the criticality value. The criticality value is the negative optimal value of~\eqref{eqn:itr_tangential}, and it is thus dependent on the norm, which we are allowed to choose in a way best suited for the problem geometry or the inner solver(s). At this point, the importance of~\cref{ass:equivalent_norms} becomes clear. It allows us to relate criticality values for different norms.

\begin{lemm}\label{thm:critical_values_equivalent}
Suppose \cref{ass:covering,ass:lipschitz_gradients_true,ass:equivalent_norms,ass:fully_linear_models} hold. Denote by $
\omegak := \omega(\vxNk; \mfk, \linFeasApproxk, \normIt{\bullet}) $ the optimal value of~\eqref{eqn:itr_tangential} and by $\omegak_2$ the optimal value if the 2-norm is used, $\omegak_2 = \omega(\vxNk; \mfk, \linFeasApproxk, \normTwo{\bullet})$. There is a constant $\CCcrit \ge 1$ such that for any $\kk \in \mathbb{N}_0$ for which the normal step exists, it holds that
\begin{equation}
\frac{1}{\CCcrit}~
\omegak
\le
\omegak_2
\le
\CCcrit~\omegak.\label{eqn:critical_values_equivalent}
\end{equation}
\end{lemm}

The proof of~\cref{thm:critical_values_equivalent} can be found in the appendix. The lemma shows that we can relate the iteration dependent inexact critical values $\omegak$ and $\omegak_2$. Furthermore, we show in the appendix that it is then sensible to assume a sufficient decrease condition for the objective surrogate functions as per~\cref{ass:sufficient_decrease}, as long as an additional (standard) assumption holds:

\begin{assu}\label{ass:hess_norm_bounded}
The norm of all model Hessians of the objective function surrogates is uniformly bounded above, i.e., there is a positive constant $\CChessF > 0$ such that for all $\kk \in \mathbb{N}_0$
\[
\norm{ \hess\mfLk_\indf (\ve {\xi})}_2 \le \CChessF
\text{ for all $\indf=1,\ldots,\fDim$, and all $\ve {\xi}\in \linFeasApproxk(\vxk) \cap \trk$.}
\]
\end{assu}

We know how to find suitable model functions that satisfy~\cref{ass:hess_norm_bounded}, including Lagrange interpolation polynomials or \acp{rbf}, so that finally we can state the sufficient decrease assumption:

\begin{assu}[Sufficient Decrease]\label{ass:sufficient_decrease}
Suppose
\cref{ass:covering,ass:lipschitz_gradients_true,ass:fully_linear_models} hold and that $\Deltak \in (0,\DeltaMax)$. Let $\scalk := \scal{\mfk}$. If $\vnk$ is compatible, $\vdk$ is a minimizer of~\eqref{eqn:itr_tangential} at $\vxk$ for $\normIt{\bullet}$, and $\omegak_2$ is the sign-flipped optimal value from~\eqref{eqn:itr_tangential} using the $2$-norm, then there is a step-length $\stepsizek \ge 0$ such that $\vxk + \vnk + \stepsizek \vdk \in \linFeasApproxk \cap \trk$ and
\begin{equation}
\scalk \args{ \vxNk} -
\scalk \args{ \vxNk + \stepsizek \vdk}
\ge
\CCsd \omegak_2
\min \set{
\frac{\omegak_2}{\CCcritDenom},
\Deltak, 1},\label{eqn:sufficient_decrease}
\end{equation}
for constants $\CCsd \in (0,1)$ and $\CCcritDenom \ge 1$.
\end{assu}

Throughout the rest of this article we use a slightly modified measure for notational convenience:

\begin{coro}[Modified Criticality Measure]\label{thm:ScaryPitchOld}
For any $\kk\in \mathbb{N}_0$ and $\omegak$ as in~\cref{thm:critical_values_equivalent}, define the criticality measure $\chik := \min \set{ 1, \omegak}$ and denote by $\chik_2$ the corresponding value, if the 2-norm is used instead of $\normIt{\bullet}$. Then $\lim_{k\to \infty} \omegak_2 = 0$ if and only if $\lim_{k\to \infty} \chik_2 = 0$, and if \cref{ass:sufficient_decrease} holds, then it also follows (with $\CCcritDenom \ge 1$) that
\begin{equation*}
\scalk \args{ \vxNk} -
\scalk\args{ \vxNk + \stepsizek \vdk}
\ge
\CCsd \chik_2
\min \set{
\frac{\chik_2}{\CCcritDenom},
\Deltak}.\label{eqn:sufficient_decrease_critMin}
\end{equation*}
\end{coro}
\section{Discussion of the Algorithm}\label{section:algorithm}

The behavior of the algorithm is determined by several additional algorithmic parameters:
\begin{center}
\begin{tabularx}{\linewidth}{cX}
\toprule Parameter(s) & Description \\
\midrule $0 < \DeltaInit \le \DeltaMax < \infty$ & initial and maximum trust-region radius
\\
$0<\gammass \le \gammas < 1 \le \gammag$ & shrinking and growing parameters for the trust-regions
\\
$\filterConst \in (0, 1)$ & filter constant in \cref{def:filter}
\\
$0 < \CCaccept \le \CCsuccess < 1$ & acceptance thresholds for trial point test
\\
$0 < \CCcritTestchi < 1$, $0 \le \CCcritTesttheta \le \CCdeltan$ & thresholds for criticality test
\\
$\kappa_\theta \in (0,1), \psi > \frac{1}{1+{\mu}}$ & threshold parameters in~\eqref{eqn:algo_model_dec}
\\
$0 < \CCBeta < \CCMu$ and $\CCalpha \in (0,1)$ & Criticality Routine threshold factors and backtracking constant
\\
$\CCsd \in (0,1)$ & \emph{s}ufficient \emph{d}ecrease constant in \cref{ass:sufficient_decrease}
\\
$\CCDelta \in (0,1), \CCmu > 0, \mu\in (0,1)$ & constants defining compatibility in \cref{def:compatible_normal_step}
\\
$\CCdeltan > 0, \, \CCubn > 0$ & existence of normal step in \cref{ass:normal_step_exists_bounded}
\\
\bottomrule
\end{tabularx}
\end{center}
\noindent{}Below, the \emph{Criticality Test} is used instead of the intuitive stopping criterion
\enquote{ $\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normIt{\bullet} } = 0$ } or \enquote{$\chik = 0$} because the surrogate models are inexact. If an iterate is nearly feasible and nearly critical for the surrogate problem, then the criticality routine is entered and the trust-region radius is reduced to make the models more precise. At a truly critical point the routine loops infinitely, but if a point is not critical for the true functions, we exit and continue with regular iterations. Our version of the routine is a bit different from the \enquote{unrolled} variant in~\cite{easonTrustRegionFilter2016}. Whilst their version is simpler, we avoid restoration in criticality loops. For further details of the Criticality Test and the Criticality Routine we refer to~\cite{connIntroductionDerivativefreeOptimization2009}.

\subsection*{Algorithm}
\begin{enumerate}
\setcounter{enumi}{-1}
\item\label[step]{algo:init} \textbf{Initialization:} Let $\kk \leftarrow 0$, $\filter \leftarrow \emptyset$ and $\vnk = \None$. Evaluate $\vf(\vxk), \vg(\vxk), \vh(\vxk)$ and compute $\thetak = \theta\args{ \vx }$. Build surrogate models $\mfk,\mgk,\mhk$ that are fully linear in $\trk$ with radius $\DeltaPrek := \Delta\itDeltaCmd{0}$ and set $\chikPre = \None$.
\item\label[step]{algo:compat_test} \textbf{Compatibility Test:}
\begin{itemize}
\item If $\vnk \ne \None$ and $\vnk$ is compatible w.r.t. $\DeltaPrek$, go to step~\ref{algo:descent}.
\item If $\vnk = \None$, try to compute $\vnk$. If $\vnk$ exists and is compatible w.r.t. $\DeltaPrek$, go to step~\ref{algo:descent}.
\end{itemize}
\item\label[step]{algo:restoration} \textbf{Restoration:} Add $\vxk$ to the filter, set $\Deltak \leftarrow \DeltaPrek, \chik \leftarrow \chikPre$, and attempt to find a \emph{restoration step} $\vrk$ and $\DeltaPreNext>0$ for which there exists a $\DeltaPre\itDeltaCmd{\kk+1}$-compatible normal step at $\vxk + \vrk$, and for which $\vxk + \vrk$ is acceptable for $\filter$. Set $\vxNext\leftarrow \vxk + \vrk$, keep $\DeltaPreNext$ and go to step~\ref{algo:updates}.
\item\label[step]{algo:descent} \textbf{Descent Step:} Compute a descent direction $\vdk$ and $\chikPre$ with~\eqref{eqn:itr_tangential} and \cref{thm:ScaryPitchOld}.
\item\label[step]{algo:crit_test} \textbf{Criticality Test:} If $\thetak < \CCcritTesttheta$
\textbf{and} $\bigl(\chikPre < \CCcritTestchi\text{ and } \DeltaPrek > \CCMu \chikPre\bigr)$, then enter the \textbf{Criticality Routine} to get $\Deltak$ and $\chik$ and update $\vdk$ and $\vnk$. Else, set $\Deltak \leftarrow \DeltaPrek$ and $\chik \leftarrow \chikPre$.
\item\label[step]{algo:acceptance_test} \textbf{Acceptance Test:} Compute a step-size $\stepsizek > 0$ such that \cref{ass:sufficient_decrease} is fulfilled. Set $\vxSk = \vxk + \vnk + \stepsizek \vdk$. Compute $\vf(\vxSk)$ and $\theta\args{ \vxSk}$.
\begin{itemize}
\item If $\vxSk$ is not acceptable for the augmented filter $\filter \cup \{ (\thetak, \scalf{\vxk}) \}$ OR
\item If
\begin{align}
\scal{\mfk}(\vxk) -
\scal{\mfk}(\vxSk) &\ge \kappa_\theta \thetak^\psi\label{eqn:algo_model_dec}
\shortintertext{AND}
\rhok :=
\frac{
\scal{\vf}(\vxk) -
\scal{\vf}(\vxSk)}{
\scal{\mfk}(\vxk) -
\scal{\mfk}(\vxSk)} &< \CCaccept,\label{eqn:algo_ratio_test}
\end{align}
\end{itemize}
keep $\vxNext \leftarrow \vxk, \vnkN \leftarrow\vnk$, choose $\DeltaPreNext\in [\gammass \Deltak, \gammas \Deltak]$, increment $\kk$ and go to step~\ref{algo:updates}.
\item\label[step]{algo:filter_test} \textbf{Filter Test:} If~\eqref{eqn:algo_model_dec} fails, include $\vxk$ in the filter ($\theta$ -iteration).
\item\label[step]{algo:iterate_update} \textbf{Iterate Updates:} Set $\vxNext \leftarrow \vxSk$ and choose
\begin{equation}
\DeltaPreNext \in
\begin{cases}
[\gammass \Deltak, \gammas \Deltak] & \text{if $\rhok < \CCsuccess$}, \\
[\Deltak, \min\{\gammag \Deltak, \DeltaMax\}] & \text{if $\rhok\ge \CCsuccess$.}
\end{cases}\label{eqn:radius_update_formula}
\end{equation}
\item\label[step]{algo:updates} \textbf{Model Updates:} Update the surrogates for $\vxNext$ and $\DeltaPreNext$. If any of the constraint function models is dependent on the trust-region radius (see~\cref{rem:radius_dependent_surrogates}), then set $\vnkN \leftarrow\None$ and $\chikNPre \leftarrow \None$. Finally, let $\kk \leftarrow \kk + 1$ and go to step~\ref{algo:compat_test}.
\end{enumerate}


Note, that in contrast to regular trust-region methods, a trial point that does not pass the test~\eqref{eqn:algo_ratio_test} may still be accepted as the next iterate, namely in a $\theta$-iteration. The restoration step $\vx\itCmd{\kk + 1} = \vxk + \vrk$ aims at reducing the constraint violation and is usually calculated by solving
\begin{equation}
\min_{\vx} \theta\args{ \vx}.
\tag{R}\label{eqn:restoration_problem}
\end{equation}
For a detailed discussion of the role of~\eqref{eqn:algo_model_dec} and the restoration procedure, we refer to~\cite{fletcherGlobalConvergenceTrustRegion2002}. We only want to note a few important properties:
\begin{itemize}
\item No feasible iterate is ever added to the filter.
\item Unless noted otherwise, we assume that the restoration procedure is always able to find a suitable restoration step.
\item There are no two successive restoration iterations.
\item Under \cref{ass:normal_step_exists_bounded}, there can only ever be a finite number of sub-iterations in the restoration procedure until a compatible normal step exists.
\end{itemize}

\subsection{Criticality Routine}\label{section:crit_routine}

The criticality routine is provided with the current models $\mfk,\mgk,\mhk$, the preliminary criticality value $\chikPre$ from algorithm step~\ref{algo:descent}, and the preliminary trust-region radius $\DeltaPrek$. All the models as well as the normal step $\vnk$ may be modified within this sub-routine, and we set $\Deltak$ and $\chik$.

When the Criticality Routine starts, we know that $\DeltaPrek > \CCMu \chikPre$ and that $\vnk$ is compatible for $\DeltaPrek$.
\begin{enumerate}
\item\label{critAlgo:init}
Let $j\leftarrow 0$ and set $\delta_j \leftarrow \DeltaPrek$, $\chi_j \leftarrow \chikPre$ and $\vn_j \leftarrow \vnk$.
\item While $\delta_j > \CCMu \chi_j$:
\begin{enumerate}
\item\label{critAlgo:radius}
Tentatively decrease the radius: $\delta_j^+ \leftarrow \CCalpha \delta_{j} = \CCalpha^{j+1} \DeltaPrek$.
\item\label{critAlgo:model_updates}
Make the models fully linear for the new radius $\delta_j^+$.
\item\label{critAlgo:compat_test}
Solve~\eqref{eqn:itr_normal} to get $\vn_j^+$ for the new models and $\delta_j^+$.\\
If $\vn_j^+$ is not compatible for $\delta_j^+$, then \texttt{BREAK}.
\item\label{critAlgo:index}
Set $j\leftarrow j+1$, $\delta_j\leftarrow \delta_{j-1}^+$ and keep $\vn_j \leftarrow \vn_{j-1}^+$.
\item\label{critAlgo:crit_update}
Solve~\eqref{eqn:itr_tangential} and use \cref{thm:ScaryPitchOld} to get the criticality measure value $\chi_j$ for the new models from step~\ref{critAlgo:model_updates} and for $\delta_j$.
\end{enumerate}
\item\label{critAlgo:finish}
Keep the updated models as $\mfk, \mgk, \mhk$, set $\chik \leftarrow \chi_j$, $\vnk \leftarrow \vn_j$ and choose
\begin{equation}
\Deltak \leftarrow \min\set{ \max\{ \delta_j, \CCBeta \chik\}, \DeltaPrek}.\label{eqn:crit_test_final_radius}
\end{equation}
\end{enumerate}

There are two ways the routine can stop. If the radius is sufficiently small compared to the criticality, i.e., $\delta_j \le \CCMu \chi_j$, or if the next prospective normal step $\vn_j^+$ is no longer compatible for the smaller radius $\delta_j^+$. The second criterion ensures (inductively) that after the routine has stopped finitely, the normal step will be compatible for $\delta_j$ and for any radius $\Delta \in [\delta_j, \DeltaMax]$. Because it always holds that $\delta_j\le \DeltaPrek$ and $\delta_j \leq \max\{ \delta_j, \CCBeta \chik\}$, it follows from~\eqref{eqn:crit_test_final_radius} that $\vnk$ is compatible for $\Deltak$.

One reason for the Criticality Test (and the Criticality Routine) is to have a lower bound on the inexact criticality if the constraint violation is sufficiently small, i.e., $\thetak \le \CCdeltan$. Then, a normal compatible normal step always exists, and the routine can stop only because of the first condition, $\delta_j \leq \CCMu \chi_j$.
\begin{lemm}\label{thm:faithfulnessPossumed}
Suppose~\cref{ass:covering,ass:lipschitz_gradients_true,ass:fully_linear_models,ass:equivalent_norms,ass:normal_step_exists_bounded} hold. Let $\kk \in \mathbb N_0$ be the iteration index. If $\thetak \le \CCdeltan$, and if the Criticality routine does not run infinitely, then
\begin{equation}
\chik \geq \min\set{\frac{\Deltak}{\CCMu}, \CCcritTestchi}.\label{eqn:primogenitsMaysters}
\end{equation}
\end{lemm}

\begin{proof}
Equation~\eqref{eqn:primogenitsMaysters} holds trivially if $\chikPre \ge \CCcritTestchi$. In this case, the Criticality Routine is not entered, we set $\chik = \chikPre$ and $\Deltak = \DeltaPrek$, and~\eqref{eqn:primogenitsMaysters} follows. Otherwise, the Criticality Routine must finish due to
\(0 < \delta_j \le \CCMu \chi_j \le \CCMu \chik\). By~\eqref{eqn:crit_test_final_radius} we have
\[
0 <
\Deltak
\le
\max\{\delta_j, \CCBeta \chik\}
\stackrel{\CCBeta < \CCMu}
\le
\max\{\delta_j, \CCMu \chik\}
\le
\CCMu \chik.
\]
Division by $\CCMu$ again implies~\eqref{eqn:primogenitsMaysters}.
\end{proof}

Finally, note that the models $\mfk, \mgk, \mhk$ are fully linear when the Criticality Routine has finished after a finite number of iterations (even if the radius is slightly increased above $\delta_j$) because of \cref{ass:fully_linear_models}.

\section{Convergence to Quasi-Stationary Points}\label{sec:quasi_convergence}

Note that the problem~\eqref{eqn:itr_tangential} defining $\chik$ is similar to the criticality problem~\eqref{eqn:DemandIdenticalOval}, but we compute the descent step starting at $\vxNk$ instead of $\vxk$, the surrogate model gradients are used to determine a model descent step within the approximated linearized feasible set and a variable vector norm $\normIt{\bullet}$ bounds the problem (which must not necessarily equal the trust-region norm $\normTrIt{\bullet}$). For the convergence analysis, the additional uncertainties are dealt with in two steps: The first part (the remainder of this section) is concerned with proving that there is a subsequence $\{\vxkl\}$ of iterates converging to a quasi-stationary point. That is, it holds that
\[
\lim_{{\ell}\to \infty} \theta(\vxkl) = 0
\quad\text{and}\quad
\lim_{{\ell}\to \infty}
\omega\args{\vxNkl; \mfkl, \linFeasApproxkl, \normInf{\bullet}} = 0.
\]
Afterwards, \cref{section:kkt_convergence} builds on this result to show convergence to actual KKT-points.

\subsection{Comparison with Related Algorithms}\label{sec:other_algos}

Most results in this section are ``translated'' from their single-objective pendants in~\cite{fletcherGlobalConvergenceTrustRegion2002} or~\cite{waltherInexactTrustregionSQPfilter2016}. The latter article is also concerned with inexact surrogates, which have error bounds that are slightly different from ours and can become exact eventually. Hence, we cannot simply apply the maximum-scalarization and be done (unfortunately). We also have to take care of the Criticality Routine and other, more subtle differences, such as the iteration-dependent norms. That is why we have decided to cite the results from~\cite{fletcherGlobalConvergenceTrustRegion2002,waltherInexactTrustregionSQPfilter2016} explicitly whenever we introduce a similar one. If we deem the proofs very similar and easily transferable, we have stated so. Otherwise, hints are given on how to adapt them or the proofs are provided wholly for the sake of completeness.

Special mention has to be made of the single-objective algorithm in~\cite{easonTrustRegionFilter2016} that is further detailed in the dissertation~\cite{easonTrustRegionFilter2018}. This algorithm also employs a Criticality Routine and uses fully linear models. By encoding the surrogate modeling of black-box components as additional constraints to the original problem and using a nonlinear (even non-quadratic) subproblem for the computation of the normal step, the convergence analysis becomes easier due to the resulting inexactness bounds. At the time of writing, the nonlinear normal step computation did not suit our particular needs, but we think it easily possible and very worthwhile to also transfer their approach to the multi-objective case.

\subsection{Final Definitions and Requirements}

To actually investigate limit points of algorithmic iteration sequences we need
\cref{ass:closed_and_bounded}, which guarantees their existence:

\begin{assu}\label{ass:closed_and_bounded}
The set $\covering$ is contained in a closed and bounded set.
\end{assu}

A detailed discussion of~\cref{ass:closed_and_bounded} and alternatives is given in~\cite{fletcherGlobalConvergenceTrustRegion2002}. If the iterates are contained in a bounded domain, the true functions (which are continuous according to~\cref{ass:lipschitz_gradients_true}) are bounded and so are their Lipschitz-continuous gradients. With~\cref{ass:fully_linear_models} the models are fully linear and have to be bounded at the iterates as well, due to $\DeltaPrek\le\DeltaMax<\infty$:
\begin{coro}\label{thm:bounded_model_gradients}
If~\cref{ass:closed_and_bounded,ass:lipschitz_gradients_true,ass:fully_linear_models} and $\DeltaPrek \le \DeltaMax < \infty$ hold, then the norm of all function and model gradients is uniformly bounded above, i.e., there is a positive constant $\CCubj > 0$ such that for
\(
\varphi \in \set{
\fL_1, \ldots, \fL_\fDim,
\gL_1, \ldots, \gL_\gDim,
\hL_1, \ldots, \hL_\hDim}
\), for all $\kk\in\mathbb{N}_0$ and
\(
\varphi\KK \in \set{
\mfLk_1, \ldots, \mfLk_\fDim,
\mgLk_1, \ldots, \mgLk_\gDim,
\mhLk_1, \ldots, \mhLk_\hDim}
\) it holds that
\begin{equation*}
\norm{\grad \varphi(\vxk)}\le \CCubj
\quad\text{and}\quad
\normInf{ \grad \varphi\KK(\vxk)} \le \CCubj.\label{eqn:def_CCubj}
\end{equation*}
\end{coro}

Throughout this section, we consider the following iteration index sets:
\begin{defi}
The set of ``successful'' iteration indices is
\(
\acceptIndices = \{ \kk: \vxNext = \vxSk \}.
\) The set of restoration indices is defined as
\(
\restoreIndices = \{ \kk:
\vnk \text{ does not exist or~\eqref{eqn:compatible_normal_step} fails, i.e., $\vnk$ is incompatible}
\}.
\) Finally, the set of filter indices (indices of iterations that modify the filter) is
\(
\filterIndices = \{
\kk: \vxk \text{ is added to the filter }
\}
\supseteq \restoreIndices.
\)
\end{defi}

Whenever we require ``the Criticality Routine to finish finitely'' in subsequent statements, then we want it to finish after a finite number of iterations and explicitly include the case that the routine is not even entered due to the Criticality Test failing in \cref{algo:crit_test} of the algorithm.

\subsection{Convergence Analysis} Similar to~\cite{waltherInexactTrustregionSQPfilter2016} we can state \enquote{accuracy requirements} that bound the linearized constraint violation of the steps $\vnk$ and $\vsk$ by the trust-region radius $\Deltak$:
\begin{lemm}[Accuracy Requirements similar to {\cite[A.2.4]{waltherInexactTrustregionSQPfilter2016}}]\label{thm:accuracy_requirements}
Suppose~\cref{ass:covering,ass:lipschitz_gradients_true,ass:fully_linear_models,ass:equivalent_norms} hold and that $\kk\in \mathbb{N}_0$ is an iteration index. Then there is a constant $\CCaccN > 0$ such that, if $\vnk$ exists as the solution to~\eqref{eqn:itr_normal}, it holds that
\begin{equation}
\max \set{
\max_\indh \abs{ \hL_\indh(\vxk) + \noDotProd{\dh_\indh(\vxk)}{\vnk} },
\max_\indg \gL_\indg(\vxk) + \noDotProd{\dg_\indg(\vxk)}{\vnk} }
\le
\CCaccN
\DeltaPrek \normTrIt{\vnk}.\label{eqn:recantsKnuts}
\end{equation}
For any $\kk \in \mathbb{N}_0$, for which $\vnk$ exists and satisfies $\normTrIt{\vnk}\le \DeltaPrek$, it also holds that
\begin{align}
\max \set{
\max_\indh
\abs{
\hL_\indh(\vxk) + \noDotProd{\dh_\indh(\vxk)}{\vnk} },
\max_\indg \gL_\indg(\vxk) + \noDotProd{\dg_\indg(\vxk)}{\vnk} } &\le \CCaccN \DeltaPrek^2,\label{eqn:accuracy_normal}
\intertext{ and if the Criticality Routine finishes finitely and $\vsk$ is the step $\vnk + \stepsizek \vdk$ with $\vdk$ computed using~\eqref{eqn:itr_tangential} and $\normTrIt{\vsk} \le \Deltak$, we have}
\max \set{
\max_\indh
\abs{
\hL_\indh(\vxk) + \noDotProd{\dh_\indh(\vxk)}{\vsk} },
\max_\indg \gL_\indg(\vxk) + \noDotProd{\dg_\indg(\vxk)}{\vsk} } &\le \CCaccS \Deltak^2.\label{eqn:accuracy_step}
\end{align}
\end{lemm}
\begin{proof}
The constraints of~\eqref{eqn:itr_normal} ensure that
\begin{equation}
\mhk(\vxk) + \dotMatProd{\DmH(\vxk)}{\vnk} = \vZ
\quad\text{and}\quad
\mgk(\vxk) + \dotMatProd{\DmG(\vxk)}{\vnk} \le \vZ.\label{eqn:centimetersHectic}
\end{equation}
Because the models are interpolating (\cref{ass:fully_linear_models}) we also have
\(\mhk(\vxk) = \vh(\vxk)\) and
\(\mgk(\vxk) = \vg(\vxk)\). Dropping the argument $\vxk$ to improve readability and then invoking the Cauchy--Schwartz inequality and the error-bounds of the models (for the preliminary radius $\DeltaPrek$), we obtain
\[
\begin{aligned}
\max_{\indh} \abs{
\hLk_\indh +
\noDotProd{{\dhk_\indh}}{\vnk} } &=
\abs{
\hLk_{\indmax{\indh}} +
\noDotProd{{\dhk_{\indmax{\indh}}}}{\vnk} - 0 }
\stackrel{\eqref{eqn:centimetersHectic}}=
\abs{
\noDotProd{
\args{
\dhk_{\indmax{\indh}} - \dmhk_{\indmax{\indh}}}}{\vnk} }
\\
&
\le
\normTwo{
\dhk_{\indmax{\indh}} - \dmhk_{\indmax{\indh}} }
\normTwo{\vnk}
\le
\CCflDiffH \DeltaPrek \CCnorm \normTrIt{\vnk},
\end{aligned}
\]
where $\indmax{\indh}$ is the maximizing index. Similarly, we find that
\[
\max_{\indg} \gL_\indg(\vxk) +
\noDotProd{\dg_\indg(\vxk)}{\vnk}
\le \CCflDiffG \DeltaPrek \CCnorm \normTrIt{\vnk},
\]
and~\eqref{eqn:recantsKnuts} follows with $\CCaccN = \CCnorm \max\{ \CCflDiffH, \CCflDiffG \}$. Moreover, $\normTrIt{\vnk} \le \DeltaPrek$ leads to~\eqref{eqn:accuracy_normal}.

The second inequality~\eqref{eqn:accuracy_step} is derived analogously, respecting the fact that the step-size is chosen so that $\normTrIt{\vsk} \le \Deltak$ and that $\vxk + \vsk$ also satisfies the approximated linearized constraints.
\end{proof}

From the preceding accuracy results, a bound on the constraint violation can be derived:

\begin{lemm}[{\cite[Lemma~4.3]{waltherInexactTrustregionSQPfilter2016}, \cite[Lemma~3.1]{fletcherGlobalConvergenceTrustRegion2002}}]\label{thm:bound_on_constraint_violation1}
Assume that the algorithm is applied to~\eqref{eqn:mop_nonlin} and that
\cref{ass:covering,ass:closed_and_bounded,ass:fully_linear_models,ass:normal_step_exists_bounded,ass:lipschitz_gradients_true,ass:equivalent_norms} hold. There is a constant $\CCubj > 0$ (independent of $\kk$) such that, if the normal step exists, it holds that
\begin{equation}
\thetak \le \args{ \CCaccN \DeltaPrek + \CCubj \CCnorm } \normTrIt{\vnk}.\label{eqn:balladedAdzukis}
\end{equation}
\end{lemm}

\begin{proof}
The proof works very similar to that of Lemma~4.3 in~\cite{waltherInexactTrustregionSQPfilter2016}. The bound (and the constant) follow from~\cref{thm:accuracy_requirements} by uniformly bounding the norm of all surrogate gradients.
\end{proof}

When the normal step is compatible, i.e., for iterations with $\kk \notin \restoreIndices$, we can further refine the bound on $\thetak$:
\begin{lemm}[{\cite[Lemma~4.4]{waltherInexactTrustregionSQPfilter2016}, \cite[Lemma~3.4]{fletcherGlobalConvergenceTrustRegion2002}}]\label{thm:bound_on_constraint_violation2}
Suppose that the algorithm is applied to~\eqref{eqn:mop_nonlin} and that
\cref{ass:covering,ass:closed_and_bounded,ass:lipschitz_gradients_true,ass:fully_linear_models,ass:normal_step_exists_bounded,ass:lipschitz_gradients_true,ass:equivalent_norms} hold. Suppose further that $\kk\notin\restoreIndices$ and that~\eqref{eqn:balladedAdzukis} holds. Then there is a constant $\CCthetaub > 0$ such that
\begin{equation}
\thetak \le \CCthetaub \max \set{ \DeltaPrek^{1+{\mu}}, \DeltaPrek^2 }
\quad\text{and}\quad
\theta\args{ \vxk + \vsk } \le \CCthetaub \Deltak^2.\label{eqn:shoalinessesDespicablenesses}
\end{equation}
\end{lemm}

\begin{proof}
Since $\kk\notin\restoreIndices$, one obtains from~\eqref{eqn:balladedAdzukis} and from~\eqref{eqn:compatible_normal_step} that
\begin{align*}
\thetak &\le
\args{ \CCaccN \DeltaPrek + \CCubj \CCnorm } \normTrIt{\vnk}
\le
\args{ \CCaccN \DeltaPrek + \CCubj \CCnorm }
\CCDelta \DeltaPrek \min \set{ 1, \CCmu \DeltaPrek^\mu }
\\
&\le
\CCDelta
\args{ \CCaccN \DeltaPrek^2 + \CCubj \CCnorm \CCmu \DeltaPrek^{1+{\mu}}}
\le
\underbrace{\CCDelta
\args{
\CCaccN + \CCubj \CCnorm \CCmu }}_{=:\CCthetaub^*\text{ (const.)}} \max \set{ \DeltaPrek^{1+\mu}, \DeltaPrek^2 }.
\end{align*}
For the second bound let $c_\ind$ be any constraint component from $\vh$ or $\vg$. Due to~\cref{ass:covering,ass:lipschitz_gradients_true} we can construct a Taylor polynomial and use the Mean-Value-Theorem to obtain
\begin{equation}
c_\ind(\vxk + \vsk) = c_\ind(\vxk) + \noDotProd{\grad c_\ind(\vxk)}{\vsk} +
\frac{
\noDotProd{{\vsk}}{
\noDotMatProd{\hess c_\ind(\ve {\xi})}{\vsk} } }{2}\label{eqn:proof_bound_on_constraint_violation2_1}
\end{equation}
for some $\ve {\xi}$ on the line segment from $\vxk$ to $\vxk + \vsk$. From~\eqref{eqn:accuracy_step} we know that
\(\abs{ c_\ind(\vxk) + \noDotProd{\grad c_\ind(\vxk)}{\vsk} }
\le
\CCaccS \Deltak^2.
\) With~\cref{ass:closed_and_bounded,ass:lipschitz_gradients_true} the norm of all constraint Hessians in
\cref{eqn:proof_bound_on_constraint_violation2_1} can be bounded from above globally (say, by $2{\kappa}_{H} > 0$) so that the Cauchy--Schwartz inequality gives
\[
\frac{\noDotProd{{\vsk}}{
\noDotMatProd{\hess c_\ind(\ve {\xi})}{\vsk} }}{2}
\le
\frac{\normTwo{\vsk}^2 \normTwo{\hess c_\ind(\ve {\xi})}}{2}
\le
\kappa_{H} \Deltak^2.
\]
Hence,~\eqref{eqn:shoalinessesDespicablenesses} follows with $
\CCthetaub =: \max \set{
\CCthetaub^*,
\CCaccS + \kappa_{H} }. $
\end{proof}

The requirements for the accuracy requirements to hold are also met within the Criticality Routine, allowing for the following corollary:
\begin{coro}
Let $\kk \in \mathbb{N}_0$ be an iteration index and suppose the same requirements as in~\cref{thm:bound_on_constraint_violation1} hold. If the Criticality Routine is entered, then for any sub-iteration index $j\in \mathbb N_0$, it holds that
\begin{equation}
\thetak
\le
\CCaccN
\delta_j^2 +
\CCubj \normTwo{\vnk_j},\label{eqn:bound_on_constraint_violation1_inloop}
\end{equation}
If the Criticality Routine is not entered or if it has completed finitely, then it also holds that
\begin{align}
\thetak &\le
\CCaccN \Deltak^2 + \CCubj \normTwo{\vnk},\label{eqn:bound_on_constraint_violation1_finCrit}
\shortintertext{and if $\kk\notin\restoreIndices$, then}
\thetak &\le
\CCthetaub \max \set{ \Deltak^{1+{\mu}}, \Deltak^2 }.\label{eqn:shoalinessesDespicablenesses_finCrit}
\end{align}
\end{coro}


\subsubsection{Convergence in the Criticality Routine}



Our first important observation is that quasi-criticality is approximated if the criticality loop runs infinitely:
\begin{lemm}[{\cite[Lemma~8]{berkemeierDerivativeFreeMultiobjectiveTrust2021}}]\label{thm:infinite_crit_routine}
Assume that
\cref{ass:covering,ass:closed_and_bounded,ass:fully_linear_models,ass:normal_step_exists_bounded,ass:lipschitz_gradients_true,ass:equivalent_norms} hold. Denote by $\chi_j$ the inexact criticality value from~\cref{thm:ScaryPitchOld} in the criticality subroutine. If the criticality routine runs infinitely ($j\to \infty$) at $\vxk$, then
\[
\theta\args{ \vxk } = 0
\text{ and }
\lim_{j\to \infty} \delta_j = 0
\text{ and }
\lim_{j\to \infty}\chi_j = 0.
\]
\end{lemm}

\begin{proof}
As in the description of the Criticality Routine, let $j \in \mathbb N_0$ be the sub-iteration index, $\delta_j^+ := \CCalpha^{j+1} \DeltaPrek$, and $\delta_j = \delta_{j}^+$ after the index $j$ has been increased. In accordance with~\cref{thm:infinite_crit_routine}, denote by $\chi_j$ the updated (doubly inexact) criticality measure from step~\ref{critAlgo:init} and step~\ref{critAlgo:crit_update}. There are two logically disjunct stopping criteria for the Criticality Routine.

Firstly, the Criticality Routine may stop if
\(
\CCMu \chi_j
\ge
\delta_j
\) for some $j\in \mathbb{N}_0$. Conversely, if the routine loops infinitely, then it must hold that
\[
\chi_j <
\frac{\delta_j}{\CCMu} =
\CCalpha^{j} \frac{\DeltaPrek}{\CCMu}
\qquad
\text{for all $j\in \mathbb{N}_0$.}
\]
Because of $\CCalpha\in (0,1)$ the right side goes to zero and then the second limit in~\cref{thm:infinite_crit_routine} also follows.

Secondly, the routine stops if no compatible step $\vn_j$ for the preliminary radius $\delta_j^+$ can be found anymore. Vice versa, if the routine runs infinitely, there always is a compatible step with
\[
\normTwo{\vn_j}
\le
\CCnorm
\CCDelta \delta_j
\min \set{ 1,
\CCmu \delta_j^\mu }.
\]
This implies $\vnk_j \to 0$ for $j\to \infty$. From~\eqref{eqn:bound_on_constraint_violation1_inloop} it then follows that it must already hold that $\theta(\vxk) = 0$.
\end{proof}

\subsubsection{Infinitely Many Filter Iterations}

Outside of the criticality loop, we can show that feasibility is approached if the number of filter iterations is infinite:

\begin{lemm}[{\cite[Lemma~3.3]{fletcherGlobalConvergenceTrustRegion2002}}]\label{thm:infinite_filter_iters_imply_feasibility}
Suppose that the algorithm is applied to~\eqref{eqn:mop_nonlin} and that
\cref{ass:covering,ass:lipschitz_gradients_true,ass:closed_and_bounded} hold. If $\abs{\filterIndices} = \infty$, then
\(
\lim_{\substack{\kk \to \infty,\\ \kk\in \filterIndices}} \thetak = 0.
\)
\end{lemm}

\begin{proof}
The proof is exactly the same as in~\cite{fletcherGlobalConvergenceTrustRegion2002}, only $f(\bullet)$ has to be substituted by $\scalf{\bullet}$. It also works with a $\fDim + 1$ dimensional filter that uses $\vf$ instead of $\scalf{\bullet}$.
\end{proof}

The next result shows that feasibility is approached for
\emph{any} subsequence of iteration indices. The result is a bit stronger than that found in~\cite{fletcherGlobalConvergenceTrustRegion2002} and not necessarily needed for the convergence proofs. A similar theorem can be found in~\cite{wachterLineSearchFilter2005}.
\begin{lemm}\label{mythm:infinite_filter_iters_imply_feasibility}
Suppose that the algorithm is applied to~\eqref{eqn:mop_nonlin} and that
\cref{ass:covering,ass:closed_and_bounded,ass:lipschitz_gradients_true} hold. If $\abs{\filterIndices} = \infty$, then
\(
\lim_{\kk \to \infty} \thetak = 0.
\)
\end{lemm}
\begin{proof}
Assume, for sake of contradiction, that there is a subsequence of indices $\{\kk_\ell\}$ for which the constraint violation is bounded below
\begin{equation}
\theta\itthetaCmd{\kk_{\ell}} \ge \varepsilon > 0\qquad \forall \kk_\ell.\label{eqn:perturbationsBibliopegy}
\end{equation}
From~\cref{thm:infinite_filter_iters_imply_feasibility} we know that there is some $\kk_0$ such that
\begin{equation}
\theta\itthetaCmd{\kk} < \varepsilon
\quad\text{for all $
\kk\in \filterIndices$ with $\kk\ge \kk_0$.}\label{eqn:humanizesTemporizers}
\end{equation}
Because there are infinitely many indices both in $\{\kk_\ell\}$ and $\filterIndices$, each $\kk\in \{\kk_\ell\}$ (except maybe the first) must lie between a smallest filter index $\kappa_1\in \filterIndices$ and a largest index $\kappa_2\in\filterIndices$, respectively, $\kappa_1(\kk) \le \kk \le \kappa_2(\kk)$, and from~\eqref{eqn:perturbationsBibliopegy} and~\eqref{eqn:humanizesTemporizers} we can deduce that there is a $\kk_1 \ge \kk_0$ such that
\[
\kk_0 < \kappa_1(\kk) < \kk < \kappa_2(\kk)
\quad \text{for all $\kk\in \{\kk_\ell\}$ with $\kk \ge \kk_1$.}
\]
Let $\kk\in \{\kk_\ell\}$ be an index with $\kk \ge \kk_1$. For all indices $\kappa$ with $\kappa_1(\kk) < \kappa < \kappa_2(\kk)$, the iteration is not a filter iteration and hence the function $\scalShort$ is ``monotonic'':
\[
\scalShort\args{\vx\itCmd{\kappa}} -
\scalShort\args{\vx\itCmd{\kappa+1}}
\ge 0.
\]
Moreover, there must be some smallest successful index $\kappa_2^*(\kk)$ with $\kk \le \kappa_2^*(\kk) < \kappa_2(\kk)$, reducing the constraint violation from above $\varepsilon$ to below $\varepsilon$ again. $\kappa_2^*(\kk)$ is not a filter index, so~\eqref{eqn:algo_model_dec} succeeds and ${\rho}\itCmd{\kappa_2^*(\kk)} \ge \CCaccept > 0$, resulting in
\begin{equation}
\scalShort\args{ \vx\itCmd{\kappa_2^*(\kk)} } -
\scalShort\args{ \vx\itCmd{\kappa_2^*(\kk) + 1} }
\ge
\CCaccept
\kappa_\theta
\theta\itthetaCmd{\kappa_2^*(\kk)}^\psi
\ge
\CCaccept
\kappa_\theta
\varepsilon^\psi > 0.\label{eqn:samovarDemiurgic}
\end{equation}
The function $\scalShort$ might be increased eventually in the filter iteration $\kappa_2(\kk)$ (when accepting $\vx\itCmd{\kappa_2(\kk) +1}$), but due to the monotonicity of the intermediate iterations it holds that
\(
\scalShort\args{\vx\itCmd{\kappa_2^*(\kk) + 1}} -
\scalShort\args{\vx\itCmd{\kappa_2(\kk)}}
\ge 0
\) and thus
\begin{equation}
\scalShort\args{ \vx\itCmd{\kappa_2^*(\kk)} } -
\scalShort\args{ \vx\itCmd{\kappa_2(\kk)} }
\ge
\CCaccept
\kappa_\theta
\varepsilon^\psi > 0,\label{eqn:bulgineOped}
\end{equation}
and $\vx\itCmd{\kappa_2(\kk)}$ is added to the filter.

For $\kk$ from above, let $\kk_+ \in \{\kk_\ell\}$ be the smallest iteration index following $\kk$ such that it holds for the enclosing filter indices that $\kappa_2(\kk) \le \kappa_1(\kk_+)$. There then must be a largest successful index $\kappa_1^*(\kk_+)$ with $\kappa_1(\kk_+) \le \kappa_1^*(\kk_+) < \kk_+$, such that the constraint violation is increased from below $\varepsilon$ to above $\varepsilon$,
\[
\theta\itthetaCmd{\kappa_1^*(\kk_+)} <
\varepsilon
\le
\theta\itthetaCmd{\kappa_1^*(\kk_+) + 1},
\]
and from the filter mechanism we know that then the objective reduction has to be significant compared to \emph{all} points in the filter with a smaller constraint violation, including that with index $\kappa_2(\kk)$ and $\theta\itthetaCmd{\kappa_2(\kk)}<\varepsilon$:
\[
\scalShort\args{ \vx\itCmd{\kappa_2(\kk)} } -
\scalShort\args{ \vx\itCmd{\kappa_1^*(\kk_+) + 1} }
\ge
\filterConst \theta\itthetaCmd{\kappa_2(\kk)}
\ge 0.
\]

Using the same notation as above, (i.e., $\kappa_2^*(\kk_+)$ as in~\eqref{eqn:samovarDemiurgic}), we deduce from the monotonicity of the intermediate iterations that
\(
\scalShort\args{ \vx\itCmd{\kappa_2(\kk)} } -
\scalShort\args{ \vx\itCmd{\kappa_2^*(\kk_+)} }
\ge 0
\) and thus, with {\eqref{eqn:bulgineOped}},
\begin{equation}
\scalShort\args{ \vx\itCmd{\kappa_2^*(\kk)} } -
\scalShort\args{ \vx\itCmd{\kappa_2^*(\kk_+)} }
\ge
\CCaccept
\kappa_\theta
\varepsilon^\psi > 0.\label{eqn:aluminousVixens}
\end{equation}

By repeating the above procedure, we see that it is possible to construct an infinite subsequence $\{\kappa_j\}$ of iteration indices from the $\kappa_2^*$ values (of indices from $\{\kk_\ell\}$), for which $\scalShort(\vx\itCmd{\kappa_j})$ is strictly monotonically decreasing with a guaranteed (constant) objective reduction~\eqref{eqn:aluminousVixens}. This is a contradiction to $\scalShort$ being bounded below as per
\cref{ass:closed_and_bounded,ass:lipschitz_gradients_true}.
\end{proof}

What follows next is a series of auxiliary lemmata to finally show convergence of the inexact criticality measure when $\abs{\filterIndices} = \infty$. First, we have a bound on the surrogate objective change along the normal step:
\begin{lemm}[{Similar to a bound in~\cite[Lemma~3.5]{fletcherGlobalConvergenceTrustRegion2002}}]\label{thm:normal_step_decrease_ub}
Suppose the algorithm is applied to~\eqref{eqn:mop_nonlin}, that
\cref{ass:covering,ass:closed_and_bounded,ass:lipschitz_gradients_true,ass:fully_linear_models,ass:equivalent_norms,ass:hess_norm_bounded} hold and that the Criticality Routine does not run infinitely. Suppose further that $\kk\notin\restoreIndices$. Then it holds that
\begin{align}
\abs{
\scalk(\vxk) - \scalk(\vxNk) } =
\abs{
\scalk(\vxNk) - \scalk(\vxk) } &\le
\CCubj \normTwo{{\vnk}} + \frac{1}{2}\CChessF \normTwo{{\vnk}}^2\label{eqn:deafenedPonytailed}
.
\end{align}
\end{lemm}

\begin{proof}
There are maximizing indices $\indf_1, \indf_2\in\{1,\ldots,\fDim\}$ such that
\[
\scalk(\vxNk) - \scalk(\vxk) =
\mfLk_{\indf_1}(\vxNk) - \mfLk_{\indf_2}(\vxk)
\le
\mfLk_{\indf_1}(\vxNk) - \mfLk_{\indf_1}(\vxk).
\]
Using a 2nd degree Taylor approximation of $\mfLk_{\indf_1}$ around $\vxk$ at $\vxNk$ results in
\begin{align*}
\scalk(\vxNk) - \scalk(\vxk) &\le
\mfLk_{\indf_1}(\vxNk) - \mfLk_{\indf_1}(\vxk)
\\
&\le
\abs{
\noDotProd{{\vnk}}{
\grad \mfLk_{\indf_1}(\vxk)} } +
\abs{
\frac{1}{2}
\noDotProd[4]{{\vnk}}{
\noDotMatProd{\hess \mfLk_{\indf_1}(\vxk)}{\vnk}} }
\\
&\le
\normTwo{\grad \mfk_{\indf_1}(\vxk)}\normTwo{\vnk} + \frac{1}{2}\normTwo{\hess \mfLk_{\indf_1}(\vxk)} \normTwo{{\vnk}}^2
\\
&
\le
\CCubj \normTwo{{\vnk}} + \frac{1}{2}\CChessF \normTwo{{\vnk}}^2,
\end{align*}
where the last inequality comes from~\cref{ass:hess_norm_bounded,thm:bounded_model_gradients}. Analogously, we can show
\begin{equation*}
-\args{\scalk(\vxNk) - \scalk(\vxk)}
\le
\mfLk_{\indf_2}(\vxk) - \mfLk_{\indf_2}(\vxNk)
\le
\CCubj \normTwo{{\vnk}} + \frac{1}{2}\CChessF \normTwo{{\vnk}}^2.\qedhere
\end{equation*}
\end{proof}

The next lemma provides a sufficient decrease bound in case that the doubly inexact criticality measure is bounded below, and the radius is sufficiently small:
\begin{lemm}[{\cite[Lemma~3.5]{fletcherGlobalConvergenceTrustRegion2002}}]\label{thm:sufficient_decrease_when_crit_bounded}
Suppose the algorithm is applied to~\eqref{eqn:mop_nonlin}, that
\cref{ass:covering,ass:closed_and_bounded,ass:lipschitz_gradients_true,ass:fully_linear_models,ass:equivalent_norms,ass:hess_norm_bounded,ass:sufficient_decrease} hold and that the Criticality Routine does not run infinitely. Suppose further, that $\kk\notin\restoreIndices$, that
\begin{equation}
\chik \ge \frac{1}{\CCcrit} \chik_2 > \frac{1}{\CCcrit} \CCepschi
\tag{LB}\label{eqn:porcellanizingTrichromat}
\end{equation}
for some $\CCepschi>0$ and that
\begin{equation}\label{eqn:pentanolFice}
\Deltak
\le
\delta_m :=
\min \Bigg\{ \;
\frac{
\CCepschi }{
\CCcritDenom }, {\args{
\frac{ 2 \CCubj }{
\CChessF\CCDelta\CCmu } }}^{\frac{1}{1+\mu}},
\args{
\frac{
\CCsd \CCepschi } { 4\CCubj \CCDelta \CCmu }^{\frac{1}{1+\mu}} }
\; \Bigg\}.
\end{equation}
Then
\begin{equation}
\scalk(\vxk) - \scalk(\vxk + \vsk)
\ge
\frac{1}{2}
\CCsd \CCepschi \Deltak.\label{eqn:foreversTalced}
\end{equation}
\end{lemm}
\begin{proof}
The proof is similar to that given in~\cite{fletcherGlobalConvergenceTrustRegion2002}. The bounds~\eqref{eqn:porcellanizingTrichromat} and~\eqref{eqn:pentanolFice} are substituted into the sufficient decrease equation of~\cref{ass:sufficient_decrease}. After some arithmetic transformations, the bound~\eqref{eqn:foreversTalced} follows from~\cref{thm:normal_step_decrease_ub},~\eqref{eqn:pentanolFice}, and the fact that $\vnk$ is compatible.
\end{proof}

Under similar conditions as in~\cref{thm:sufficient_decrease_when_crit_bounded} the iteration will be successful:
\begin{lemm}[{\cite[Lemma~3.6]{fletcherGlobalConvergenceTrustRegion2002}}]\label{thm:success_when_crit_bounded}
Suppose the algorithm is applied to~\eqref{eqn:mop_nonlin}, that
\cref{ass:covering,ass:closed_and_bounded,ass:lipschitz_gradients_true,ass:fully_linear_models,ass:equivalent_norms,ass:hess_norm_bounded,ass:sufficient_decrease} hold and that the Criticality Routine does not run infinitely. Suppose further, that $\kk\notin\restoreIndices$, that~\eqref{eqn:porcellanizingTrichromat} holds again and that
\begin{equation}
\Deltak \le \delta_{\rho} :=
\min
\set{
\delta_m,
\frac{ (1-\CCsuccess) \CCsd \CCepschi } { 2 \CCflF } }.\label{eqn:yogurt_plant_cup}
\end{equation}
Then iteration $\kk$ is successful, that is,
\(
\rhok \ge \CCsuccess.
\)
\end{lemm}
\begin{proof}
We can again follow the proof in~\cite{fletcherGlobalConvergenceTrustRegion2002}, because from~\cref{ass:fully_linear_models} we can conclude (cf.~\cite[Lemma~4.16]{thomannTrustRegionAlgorithmHeterogeneous2019}) that the model error bound holds also for the scalarization:
\begin{equation}
\abs{
\scal{\vf}(\ve {\xi}) -
\scal{\mfk}(\ve {\xi}) } =
\abs{
\scalShort(\ve {\xi}) -
\scalk(\ve {\xi}) } \le \CCflF \Deltak^2 \quad \forall \ve {\xi} \in \trk.\label{eqn:fully_linear_scalarization}
\end{equation}
\end{proof}

Further, if the radius is sufficiently small, then the test~\eqref{eqn:algo_model_dec} will succeed (prohibiting a $\theta$-iteration):
\begin{lemm}[{\cite[Lemma~3.7]{fletcherGlobalConvergenceTrustRegion2002}}]\label{thm:model_dec_test_succeeds}
Suppose the algorithm is applied to~\eqref{eqn:mop_nonlin}, that
\cref{ass:covering,ass:closed_and_bounded,ass:normal_step_exists_bounded,ass:fully_linear_models,ass:lipschitz_gradients_true,ass:equivalent_norms,ass:hess_norm_bounded,ass:sufficient_decrease} hold and that the Criticality Routine does not run infinitely. Suppose further, that $\kk\notin\restoreIndices$, that~\eqref{eqn:porcellanizingTrichromat} holds again, that $\thetak \le \CCdeltan$, and that
\begin{equation}
\Deltak \le
\delta_f :=
\min
\set{
\delta_m, 1, {\args{
\frac{
\CCsd \CCepschi }{ 2{\kappa}_\theta \CCthetaub^{{\psi}} } }}^{\frac{1}{{\psi}(1 + \mu)-1}} }.\label{eqn:edison_tesla_boxfight}
\end{equation}
Then~\eqref{eqn:algo_model_dec} succeeds, i.e.,
\(
\scalk(\vxk) -
\scalk(\vxSk)
\ge
\kappa_\theta \thetak^\psi.
\)
\end{lemm}
\begin{proof}
\cref{thm:sufficient_decrease_when_crit_bounded} provides
\begin{equation*}
\scalk(\vxk) -
\scalk(\vxSk) =
\scalk(\vxk) - \scalk(\vxk + \vsk)
\ge
\frac{1}{2}
\CCsd \CCepschi \Deltak.
\end{equation*}
The small constraint violation $\thetak \le \CCdeltan$ and $\kk\notin \restoreIndices$ imply~\eqref{eqn:shoalinessesDespicablenesses_finCrit}. Because of $\Deltak \le 1$ and $\mu \in (0,1)$ the bound~\eqref{eqn:shoalinessesDespicablenesses_finCrit} simplifies to $\thetak \le \Deltak^{1+{\mu}}$, and it follows that
\(
\CCthetaub^{{\psi}} \Deltak^{{\psi}(1+\mu)}
\ge \thetak^\psi.
\) The final inequality in~\eqref{eqn:edison_tesla_boxfight} then leads to
\[
\frac{1}{2}
\CCsd \CCepschi \Deltak
\ge
\kappa_\theta \CCthetaub \Deltak^{{\psi}(1+\mu)},
\]
which then implies~\eqref{eqn:algo_model_dec}.
\end{proof}

If additionally the constraint violation is small enough, then the filter pair at the trial point will not be filter-dominated by the pair at $\vxk$:
\begin{lemm}[{\cite[Lemma~3.8]{fletcherGlobalConvergenceTrustRegion2002}}]\label{thm:small_radius_small_constraint_ensure_filter_acceptance}
Suppose the algorithm is applied to~\eqref{eqn:mop_nonlin}, that
\cref{ass:covering,ass:closed_and_bounded,ass:normal_step_exists_bounded,ass:fully_linear_models,ass:lipschitz_gradients_true,ass:equivalent_norms,ass:hess_norm_bounded,ass:sufficient_decrease} hold and that the Criticality Routine does not run infinitely. Suppose further, that $\kk\notin\restoreIndices$, that~\eqref{eqn:porcellanizingTrichromat} holds again, that $\Deltak \le \delta_\rho$ as in~\eqref{eqn:yogurt_plant_cup} and that $\thetak \le \CCdeltan$ and that
\begin{equation}
\thetak
\le
\delta_\theta :=
\min \set{
\frac{1}{\CCthetaub} {\args{
\frac{
\CCsuccess\CCsd\CCepschi }{ 2 \filterConst } }}^{2},
\frac{1}{\CCthetaub^{\frac{1}{\mu}}} {\args{
\frac{
\CCsuccess\CCsd\CCepschi }{ 2 \filterConst } }}^{\frac{1 + \mu}{\mu}} }.\label{eqn:fan_stapler_palm}
\end{equation}
Then
\(
\scalShort\args{\vxk} -
\scalShort\args{\vxk + \vsk}
\ge
\filterConst \thetak.
\)
\end{lemm}
\begin{proof}
From~\cref{thm:bound_on_constraint_violation2} we deduce that
\(
\frac{
\thetak }{
\CCthetaub }
\le
\Deltak^{1+\mu}
\) or
\(
\frac{
\thetak }{
\CCthetaub }
\le \Deltak^2
\) and thus
\begin{equation}
\Deltak
\ge
\min\set{ {\args{
\frac{
\thetak }{
\CCthetaub } }}^{\frac{1}{1+\mu}}, {\args{
\frac{
\thetak }{
\CCthetaub}}}^{\frac{1}{2}}}.\label{eqn:fan_stapler_palm1}
\end{equation}
With~\cref{thm:success_when_crit_bounded}, it follows that
\begin{align*}
\scalShort\args{\vxk} -
\scalShort\args{\vxk + \vsk} &
\ge
\CCsuccess \args{
\scalk\args{\vxk} -
\scalk\args{\vxk + \vsk}}
\stackrel{\mathclap{\eqref{eqn:foreversTalced}}}
\ge
\frac{1}{2}
\CCsuccess
\CCsd \CCepschi \Deltak
\\
&
\stackrel{\mathclap{\eqref{eqn:fan_stapler_palm1}}}
\ge
\frac{1}{2}
\CCsuccess
\CCsd \CCepschi
\min\set{ {\args{
\frac{
\thetak}{
\CCthetaub}}}^{\frac{1}{1+\mu}}, {\args{
\frac{
\thetak}{
\CCthetaub}}}^{\frac{1}{2}}}
\stackrel{\mathclap{\eqref{eqn:fan_stapler_palm}}}
\ge
\filterConst \thetak.\qedhere
\end{align*}
\end{proof}

In the preceding lemmata, it has always been assumed that the iteration index does not belong to the set of restoration indices $\restoreIndices$. This is ensured if both the radius and the constraint violation are small enough:
\begin{lemm}[{\cite[Lemma~3.9]{fletcherGlobalConvergenceTrustRegion2002}}]\label{thm:non_restore_iteration}
Suppose the algorithm is applied to~\eqref{eqn:mop_nonlin}, that
\cref{ass:covering,ass:closed_and_bounded,ass:normal_step_exists_bounded,ass:fully_linear_models,ass:lipschitz_gradients_true,ass:equivalent_norms,ass:hess_norm_bounded,ass:sufficient_decrease} hold and that the Criticality Routine does not run infinitely. Suppose that~\eqref{eqn:porcellanizingTrichromat} holds again. Suppose further that
\begin{equation}\label{eqn:GrabWoodenEvent}
\Deltak \le
\delta_{\restoreIndices} :=
\min \Bigg\{
\gammass \delta_\rho, {\args{
\frac{1}{
\CCmu}}}^{\frac{1}{\mu}}, {\args{
\frac{ (1 - \filterConst)\gammass^2\CCDelta\CCmu}{
\CCubn \CCthetaub}}}^{\frac{1}{1-\mu}},
\Bigg\}
\end{equation}
and that
\begin{equation}
\thetak \le \min \set{ \delta_\theta,
\CCdeltan}.\label{eqn:stoopeAquadromes}
\end{equation}
If $\kk > 0$, then $\kk \notin \restoreIndices$.
\end{lemm}

\begin{proof}\ 
The proof works exactly as in~\cite{fletcherGlobalConvergenceTrustRegion2002}, because the Criticality Routine is not entered in restoration\linebreak iterations.
\end{proof}

The previous lemmata allow us to investigate two mutually exclusive cases defined by the number of filter iterations. The first convergence result is obtained for the case that there are infinitely many such iterations.

\begin{lemm}\label{mythm:del_to_zero}
Suppose the algorithm is applied to~\eqref{eqn:mop_nonlin}, that
\cref{ass:covering,ass:closed_and_bounded,ass:normal_step_exists_bounded,ass:fully_linear_models,ass:lipschitz_gradients_true,ass:equivalent_norms,ass:hess_norm_bounded,ass:sufficient_decrease} hold and that the Criticality Routine does not run infinitely. Suppose $\abs{\filterIndices} = \infty$. For any subsequence of iteration indices $\{\kk_\ell\}$ with $\abs{ \{ \kk_\ell \} \cap \filterIndices} = \infty$ it holds that
\(
\liminf_{{\ell} \to \infty} \Deltakl = 0.
\) Especially, we have
\(
\lim_{\substack{\kk \to \infty, \\ \kk \in \filterIndices}}
\Deltak = 0.
\)
\end{lemm}


\begin{proof}
The proof is along the lines of~\cite[Lemma~3.10]{fletcherGlobalConvergenceTrustRegion2002}, but additionally takes into account the Criticality Routine. First, let $\{\kk_\ell\}$ a sequence of indices containing infinitely many filter indices, $\abs{\{\kk_\ell \}\cap \filterIndices} = \infty$, and assume that it was bounded away from zero:
\begin{equation}
\Deltakl
\ge
\DeltaMin > 0
\qquad \forall{\ell}\in \mathbb{N}_0.\label{eqn:outthrowGonadectomy}
\end{equation}
From
\cref{mythm:infinite_filter_iters_imply_feasibility} and
\cref{ass:normal_step_exists_bounded} it follows that
\begin{equation}
\lim_{{\ell}\to \infty} \normTrIt{ \vnkl} = 0
\stackrel{\text{\cref{ass:equivalent_norms}}}=
\lim_{{\ell}\to \infty} \normTwo{ \vnkl}.\label{eqn:siltsWeighted}
\end{equation}
Consequently, for sufficiently large $\ell$, the steps will become compatible and satisfy~\eqref{eqn:compatible_normal_step} and thus $\kk_\ell\notin\restoreIndices$. Hence, \cref{thm:normal_step_decrease_ub} applies and inequality~\eqref{eqn:deafenedPonytailed} holds again:
\[
\abs{
\scalkl\args{\vxNkl} -
\scalkl\args{\vxkl}}
\le \CCubj \normTwo{{\vnkl}} + \frac{1}{2}\CChessF
\normTwo{{\vnkl}}^2,
\]
so that~\eqref{eqn:siltsWeighted} gives
\begin{equation}
\lim_{{\ell}\to \infty}
\abs{
\scalkl\args{ \vxkl} -
\scalkl\args{ \vxNkl}} = 0.\label{eqn:acceptedlyUnprofessed}
\end{equation}
\cref{mythm:infinite_filter_iters_imply_feasibility} enables the Criticality Routine. Additionally, from the boundedness of $\Deltakl$ as per~\eqref{eqn:outthrowGonadectomy} and from~\cref{ass:normal_step_exists_bounded}, it also follows that, for large $\ell$, the Criticality Routine will not exit because no compatible normal step exists anymore. This ensures that
\begin{equation}
\chikl
\ge \min \set{ \CCcritTestchi, \frac{\Deltakl}{\CCMu}}
\ge \min \set{ \CCcritTestchi, \frac{\DeltaMin}{\CCMu}} =:
\mathtt{z} > 0\label{eqn:limacineAgrestal}
\end{equation}
for large $\ell$. We can use this fact in
\cref{ass:sufficient_decrease} and obtain
\begin{equation}
\scalkl \args{\vxNkl} -
\scalkl\args{\vxSkl}
\ge
\underbrace{
\CCsd
\mathtt z
\min \set{
\frac{\mathtt z}{\CCcritDenom},
\DeltaMin}}_{\text{const.}} > 0.\label{eqn:restrictivenessPassage}
\end{equation}
If we \enquote{add zero}, then we obtain the decomposition
%~\eqref{eqn:curculiosToutie}
\begin{equation}
\scalkl\args{\vxkl} -
\scalkl\args{\vxSkl} =
\args{
\scalkl\args{\vxkl} -
\scalkl\args{\vxNkl}} +
\args{
\scalkl\args{\vxNkl} -
\scalkl\args{\vxSkl}}.\label{eqn:curculiosToutie}
\end{equation}
By looking at~\eqref{eqn:acceptedlyUnprofessed}, we see that the first set of parenthesis vanishes, and that difference of values is the same in the limit:
\begin{equation}\label{eqn:synedriaBuckshish}
\lim_{{\ell}\to \infty}
\scalkl\args{\vxkl} -
\scalkl\args{\vxSkl} =
\lim_{{\ell}\to \infty}
\scalkl\args{\vxNkl} -
\scalkl\args{\vxSkl}.
\end{equation}
Because of $\abs{\{\kk_\ell\}\cap\filterIndices} = \infty$, there is a subsequence $\{\kk_j\}\subseteq \{\kk_\ell\}$ of filter indices, $\kk_j\in \filterIndices$, for which it must hold that $\kk_j \in \restoreIndices$ or that~\eqref{eqn:algo_model_dec} fails. We have already shown that $\kk_j\notin\restoreIndices$ for large $j$. Hence, for $j$ sufficiently large, it must hold that
\(
\kappa_\theta \theta\itthetaCmd{\kk_j}^\psi >
\scalShort\itCmd{\kk_j}(\vx\itCmd{\kk_j}) -
\scalShort\itCmd{\kk_j}(\vx\itCmd{\kk_j} + \vs\itCmd{\kk_j})
\) and~\cref{mythm:infinite_filter_iters_imply_feasibility} implies that the \ac{lhs} must go to zero. So the limits in~\eqref{eqn:synedriaBuckshish} are bounded above by $0$. But this is a contradiction to~\eqref{eqn:restrictivenessPassage}. Hence, no infinite subsequence $\{\kk_\ell\}$ with $\abs{\{\kk_\ell\}\cap \filterIndices} = \infty$ and~\eqref{eqn:outthrowGonadectomy} can exist.
\end{proof}

\begin{rema}\label{remark:durnedestStonen}
Without further assumptions on the restoration step it does not seem possible to show $\lim_{\kk\to \infty} \Deltak = 0$ in the case that $\abs{\filterIndices} = \infty$. For any element of an index sequence $\{\kk_\ell\}$, bounded away from zero and with $\abs{\{\kk_\ell\}\cap\filterIndices} < \infty$, there might be a preceding restoration iteration with arbitrarily small radius itself, but possibly increasing it without any restriction, so that no contradiction can be derived. This difficulty has been observed in the literature before, see e.g.~\cite[Remark~7]{wachterLineSearchFilter2005}.
\end{rema}

\begin{lemm}\label{mythm:del_bounded}
Suppose the algorithm is applied to~\eqref{eqn:mop_nonlin}, that
\cref{ass:covering,ass:closed_and_bounded,ass:normal_step_exists_bounded,ass:fully_linear_models,ass:lipschitz_gradients_true,ass:equivalent_norms,ass:hess_norm_bounded,ass:sufficient_decrease} hold and that the Criticality Routine does not run infinitely. Suppose that $\abs{\filterIndices} = \infty$ and that the doubly inexact criticality is bounded as per~\eqref{eqn:porcellanizingTrichromat} for all $\kk \in \filterIndices$. Then the trust-region radius is also bounded below, i.e., there is some $\DeltaMin > 0$ such that
\(
\Deltak \ge \DeltaMin
\) for all $k\in \filterIndices$.
\end{lemm}

\begin{proof}
If there are infinitely many filter indices $\filterIndices$, then we know from~\cref{mythm:infinite_filter_iters_imply_feasibility} that $\thetak \to 0$. Suppose ${\{\Deltak\}}_{\kk \in \filterIndices}$ is not bounded away from zero. Then there is a subsequence $\{\kk_\ell\}\subseteq \filterIndices$ with $\lim_{{\ell}\to \infty} \Deltakl = 0$. Following the argumentation in~\cite[Lemma~3.10]{fletcherGlobalConvergenceTrustRegion2002}, we see that for $\ell$ large enough,
\cref{thm:non_restore_iteration} applies and guarantees that $\kk_\ell\notin\restoreIndices$. At the same time, \cref{thm:model_dec_test_succeeds} applies for large $\ell$, so that the test~\eqref{eqn:algo_model_dec} always succeeds for $\vxkl$. Hence, there is some $\ell_0\in \mathbb{N}_0$ such that for all $\kk_\ell$ with $\ell\ge \ell_0$, it holds that $\kk_\ell \notin \filterIndices$. This contradicts $\{\kk_\ell\}\subseteq \filterIndices$.
\end{proof}

We are finally able to state the convergence result for the case of infinitely many filter iterations as a cumulative corollary derived from~\cref{thm:infinite_filter_iters_imply_feasibility,mythm:del_to_zero,mythm:del_bounded}:

\begin{coro}[{\cite[Lemma~3.10]{fletcherGlobalConvergenceTrustRegion2002}}]\label{thm:infinite_filter_iterations_convergence}
Suppose the algorithm is applied to~\eqref{eqn:mop_nonlin}, that
\cref{ass:covering,ass:closed_and_bounded,ass:normal_step_exists_bounded,ass:fully_linear_models,ass:lipschitz_gradients_true,ass:equivalent_norms,ass:hess_norm_bounded,ass:sufficient_decrease} hold and that the Criticality Routine does not run infinitely. Suppose further, that $\abs{\filterIndices} = \infty$. Then there is a subsequence $\{\kk_\ell\} \subseteq \filterIndices$ of filter indices with
\begin{equation*}\label{eqn:backpack_mug_phone1}
\lim_{{\ell}\to \infty} \theta\itthetaCmd{\kk_{\ell}} = 0,
\quad
\lim_{\ell \to \infty} \Deltakl = 0
\quad \text{and} \quad
\lim_{{\ell}\to \infty} \chi\itCmd{\kk_{\ell}} = 0.
\end{equation*}
\end{coro}

\subsubsection{Finitely Many Filter Iterations}

We have shown that a quasi-stationary point is approached when there are infinitely many filter indices. We now concentrate on the case that there are only finitely many filter indices (and note that this then implies that there are only finitely many restoration indices as well).
\emph{From now on, $\kk_0$ is the last iteration index for which $\vx\itCmd{\kk_0 -1}$ is added to the filter.}

\begin{lemm}[{\cite[Lemma~3.11]{fletcherGlobalConvergenceTrustRegion2002}}]\label{thm:finite_filter_iter_implies_feasibility}
Suppose the algorithm is applied to~\eqref{eqn:mop_nonlin}, that
\cref{ass:covering,ass:closed_and_bounded,ass:normal_step_exists_bounded,ass:fully_linear_models,ass:lipschitz_gradients_true,ass:equivalent_norms,ass:hess_norm_bounded,ass:sufficient_decrease} hold and that the Criticality Routine does not run infinitely. Suppose further, that $\abs{\filterIndices} < \infty$. Then
\begin{equation}
\lim_{\kk \to \infty} \thetak = 0.\label{eqn:delicious_gummy_bears0}
\end{equation}
Furthermore, $\vnk$ satisfies~\eqref{eqn:balladedAdzukis} for all $\kk\ge \kk_0$ large enough.
\end{lemm}

\begin{proof}
First consider the case that there are only finitely many successful indices $\acceptIndices$. Then $\Deltak \xrightarrow{\kk\to \infty} 0$. Because there are no restoration iterations for $\kk \ge \kk_0$, it follows from~\eqref{eqn:compatible_normal_step} that $\vnk \to 0$. Equation~\eqref{eqn:bound_on_constraint_violation1_finCrit} finally gives $\thetak \to 0$.

Now, suppose that $\abs{\acceptIndices} = \infty$ and consider any \emph{successful} iteration with $\kk\ge\kk_0$. Then $\vxk$ is not added to the filter, and it follows from step~\ref{algo:acceptance_test} of the algorithm that $\rhok \ge \CCaccept$ and from the definition of $\rhok$ and from~\eqref{eqn:algo_model_dec} that
\begin{equation}
\scalShort\args{\vxk} -
\scalShort\args{\vxk + \vsk}
\ge
\CCaccept
\args{
\scalk\args{\vxk} -
\scalk\args{\vxk + \vsk}}
\ge
\CCaccept \kappa_\theta \thetak^\psi \ge 0.\label{eqn:delicious_gummy_bears}
\end{equation}
The sequence ${\set{\scalShort\args{\vxk}}}_{\kk \ge \kk_0}$ is bounded below due to~\cref{ass:closed_and_bounded,ass:lipschitz_gradients_true}, and it is monotonically decreasing because~\eqref{eqn:algo_model_dec} always succeeds due to $\kk \ge \kk_0$. Thus, no new iterate is ever chosen with $\rhok < \CCaccept$. Hence, for the \ac{rhs} in~\eqref{eqn:delicious_gummy_bears} it follows that
\begin{equation}
\lim_{\substack{\kk\in\acceptIndices, \\ \kk\to \infty}}
\scalShort\args{\vxk} -
\scalShort\args{\vxk + \vsk} = 0.\label{eqn:delicious_gummy_bears1}
\end{equation}
The limit~\eqref{eqn:delicious_gummy_bears0} follows from the \ac{rhs} in~\eqref{eqn:delicious_gummy_bears} by noticing that $\theta\iterSubCmd{j} = \thetak$ for all non-successful indices $j\ge k$ with $\rho\itCmd{j}<\CCaccept$. The bound~\eqref{eqn:balladedAdzukis} holds eventually because $\thetak \le \CCdeltan$ for large $\kk$ and then~\cref{thm:bound_on_constraint_violation1} applies.
\end{proof}

The next auxiliary lemma shows that the trust-region radius is bounded below if the asymptotically feasible iterates do not approach a quasi-stationary point. It is used afterwards to derive a contradiction.
\begin{lemm}[{\cite[Lemma~3.12]{fletcherGlobalConvergenceTrustRegion2002}}]\label{thm:finite_filter_iterations_bounded_radius}
Suppose the algorithm is applied to~\eqref{eqn:mop_nonlin}, that
\cref{ass:covering,ass:closed_and_bounded,ass:normal_step_exists_bounded,ass:fully_linear_models,ass:lipschitz_gradients_true,ass:equivalent_norms,ass:hess_norm_bounded,ass:sufficient_decrease} hold and that the Criticality Routine does not run infinitely. Suppose further, that $\abs{\filterIndices} < \infty$ and that the doubly inexact criticality is bounded away from zero as in~\eqref{eqn:porcellanizingTrichromat} for all $\kk\ge\kk_0$. Then there is a $\DeltaMin > 0$ such that
\(
\Deltak \ge \DeltaMin
\) for all $k \in \mathbb{N}_0$.
\end{lemm}
\begin{proof}
We sketch how to adapt the proof in~\cite{fletcherGlobalConvergenceTrustRegion2002}:
\cref{thm:finite_filter_iter_implies_feasibility} ensures that for large enough $\kk$ we have $\thetak \le \CCcritTesttheta$ and the Criticality Routine affects the criticality value and the radius. For these $\kk$ it follows from~\eqref{eqn:crit_test_final_radius} and~\eqref{eqn:porcellanizingTrichromat} that
\(
\Deltak
\ge
\min\set{ \CCBeta \omegak, \DeltaPrek}
\ge
\min\set{ \CCBeta \CCepschi, \DeltaPrek}
\). Further, if $\kk \ge 1$ and $\kk -1\notin\restoreIndices$, then
\(
\Deltak
\ge
\min\set{ \CCBeta \CCepschi, \gammass \Delta\itDeltaCmd{\kk-1}}.
\) Suppose $\kk_1\ge\kk_0$ is large enough that $\thetak \le \CCcritTesttheta$ and~\eqref{eqn:stoopeAquadromes} is fulfilled for all $\kk\ge \kk_1$, i.e., $\thetak \le \min\{ \CCcritTesttheta, \delta_\theta \}$ for all $\kk\ge \kk_1$. For the purpose of deriving a contradiction, one may now assume that $j \ge \kk_1$ is the first index with
\begin{equation}
\Delta\itDeltaCmd{j}
\le
\gammass \min \set{
\delta_\rho,
\sqrt{
\frac{ (1- \filterConst)\theta^{\filter}}{
\CCthetaub}},
\Delta\itDeltaCmd{\kk_1},
\CCBeta \CCepschi} =:
\gammass \delta_s,\label{eqn:nicenessesThreadiest}
\end{equation}
where $\delta_\rho$ is as defined in~\cref{thm:success_when_crit_bounded} and
\(
\theta^{\filter} :=
\min_{i\in \filterIndices} \theta\iterSubCmd{i}
\) and proceed as in~\cite{fletcherGlobalConvergenceTrustRegion2002}.
\end{proof}

\begin{lemm}[{\cite[Lemma~3.13]{fletcherGlobalConvergenceTrustRegion2002}}]\label{mythm:finite_filter_iterations_radius_to_zero}
Suppose the algorithm is applied to~\eqref{eqn:mop_nonlin}, that
\cref{ass:covering,ass:closed_and_bounded,ass:normal_step_exists_bounded,ass:fully_linear_models,ass:lipschitz_gradients_true,ass:equivalent_norms,ass:hess_norm_bounded,ass:sufficient_decrease} hold and that the Criticality Routine does not run infinitely. Suppose further, that $\abs{\filterIndices} < \infty$. Then
\(
\liminf_{\kk \to \infty} \Deltak = 0.
\)
\end{lemm}

\begin{proof}
If there are only finitely many successful indices, then the result follows immediately from the radius update rules.

Thus, suppose that $\abs{\acceptIndices} = \infty$ and, for sake of contradiction, assume that there is a constant $\DeltaMin > 0$ such that
\begin{equation}
\Deltak \ge \DeltaMin > 0 \qquad \forall\kk\in \mathbb{N}_0.\label{eqn:jejunityThereanent}
\end{equation}
As before, we see that from~\cref{thm:finite_filter_iter_implies_feasibility} it follows for large $\kk$ that the normal step $\vnk$ exists and satisfies~\eqref{eqn:normal_step_exists_bounded} in~\cref{ass:normal_step_exists_bounded}. Then the same equations as in~\cref{thm:finite_filter_iter_implies_feasibility} hold, namely~\eqref{eqn:delicious_gummy_bears} and~\eqref{eqn:delicious_gummy_bears1}. Very much like in the proof of~\cref{thm:infinite_filter_iterations_convergence}, we can again decompose the model decrease via~\eqref{eqn:curculiosToutie}. For the first term, \eqref{eqn:deafenedPonytailed} applies again, yielding
\(
\lim_{\kk \to \infty}
\args{
\scalk\args{\vxk} -
\scalk\args{\vxNk}} = 0.
\) From equations~\eqref{eqn:delicious_gummy_bears},~\eqref{eqn:delicious_gummy_bears1} and from the model decrease decomposition it follows that
\begin{equation}
\lim_{
\substack{
\kk \to \infty, \\
\kk \in \acceptIndices}}
\args{
\scalk\args{\vxNk} -
\scalk\args{\vxSk}} =0.\label{eqn:slowingsRecreance}
\end{equation}
But the sufficient decrease condition (\cref{ass:sufficient_decrease}) still holds and (for large $\kk$ and assuming~\eqref{eqn:jejunityThereanent} holds) the Criticality Routine ensures
\[
\chik
\ge
\min\set{ \CCcritTestchi, \frac{\Deltak}{\CCMu}}
\ge
\min\set{ \CCcritTestchi, \frac{\DeltaMin}{\CCMu}} =: \mathtt{z} > 0,
\]
which can be plugged in and gives:
\begin{align*}
\scalk\args{\vxNk} -
\scalk\args{\vxSk}
\ge
\CCsd
\mathtt{z}
\min \set{
\frac{\mathtt{z}}{\CCcritDenom},
\DeltaMin} > 0,
\end{align*}
where the \ac{rhs} is constant, contradicting~\eqref{eqn:slowingsRecreance}.
\end{proof}

\begin{coro}[{\cite[Lemma~3.13]{fletcherGlobalConvergenceTrustRegion2002}}]\label{thm:infinite_filter_iterations_crit_min_convergence}
Suppose the algorithm is applied to~\eqref{eqn:mop_nonlin}, that
\cref{ass:covering,ass:closed_and_bounded,ass:normal_step_exists_bounded,ass:fully_linear_models,ass:lipschitz_gradients_true,ass:equivalent_norms,ass:hess_norm_bounded,ass:sufficient_decrease} hold and that the Criticality Routine does not run infinitely. Suppose further, that $\abs{\filterIndices} < \infty$. Then
\[
\liminf_{\kk\to \infty} \chik = 0.
\]
For a subsequence $\{\kk_\ell\} \subseteq \{\kk\}$ with $\chikl \to 0$ it must also hold that $\lim_{{\ell}\to \infty} \Deltakl = 0$.
\end{coro}

\begin{proof}
To derive a contradiction, it is assumed that $\chik$ is bounded away from zero as per~\eqref{eqn:porcellanizingTrichromat} for all $\kk \in \mathbb{N}_0$. Then the radius must be bounded away from zero due to
\cref{thm:finite_filter_iterations_bounded_radius}, in contradiction to $\liminf_{\kk\to \infty} \Deltak= 0$, as assured by
\cref{mythm:finite_filter_iterations_radius_to_zero}. If $\{\kk_\ell\}$ is such that $\chikl \to 0$, but we assume $\lim_{{\ell}\to \infty} \Deltakl = \DeltaMin > 0$, then the Criticality Routine again ensures~\eqref{eqn:limacineAgrestal} for large $\ell$, which is a contradiction.
\end{proof}

\begin{theo}[Convergence to Quasi-Stationary Point]\label{thm:convergence1}
Suppose the algorithm is applied to~\eqref{eqn:mop_nonlin} and that
\cref{ass:covering,ass:closed_and_bounded,ass:normal_step_exists_bounded,ass:fully_linear_models,ass:lipschitz_gradients_true,ass:equivalent_norms,ass:hess_norm_bounded,ass:sufficient_decrease} hold. Let $\{\vxk\}$ be the sequence of iterates produced by the algorithm. Then either the restoration procedure in step~\ref{algo:restoration} terminates unsuccessfully by converging to an infeasible, first-order critical point of problem~\eqref{eqn:restoration_problem} or there is a subsequence $\{\kk_\ell\}$ of indices for which
\(
\lim_{{\ell}\to \infty} \vxkl = \bar{\vx}
\) and $\bar{\vx}$ is quasi-stationary, i.e.,
\[
\theta\args{ \bar{\vx}} =
\lim_{{\ell}\to \infty} \thetakl = 0
\text{ and }
\lim_{{\ell}\to \infty} \chikl = 0,
\]
and for which it also holds that $\lim_{{\ell}\to \infty} \Deltakl = 0$.
\end{theo}
\begin{proof}
If the restoration procedure always terminates successfully, then the convergence to a quasi-stationary point follows from~\cref{ass:closed_and_bounded} and
\cref{thm:infinite_filter_iterations_convergence,thm:infinite_crit_routine,thm:infinite_filter_iters_imply_feasibility,thm:infinite_filter_iterations_crit_min_convergence}.
\end{proof}



\section{Convergence to KKT-Points}\label{section:kkt_convergence}

In this section, we conclude that a quasi-stationary point is also a KKT-point under suitable constraint qualifications. In the first step, we show that we can easily replace the
\emph{objective} surrogates with their true counterparts.
\begin{lemm}\label{thm:true_objective_criticality}
Suppose \cref{ass:covering,ass:lipschitz_gradients_true,ass:fully_linear_models,ass:normal_step_exists_bounded,ass:equivalent_norms} hold and suppose that it holds for a subsequence $\{\vxk\}$ of iterates that
\[
\lim_{\kk \to \infty}\Deltak = 0
\quad\text{and}\quad
\lim_{\kk \to \infty} \omegak_2 =
\lim_{\kk \to \infty} \omega\args{ \vxNk; \mfk, \linFeasApproxk, \normTwo{\bullet}} = 0.
\]
Then it also holds for the true objectives $\vf$ that
\[
\lim_{\kk \to \infty} \omegak_2 =
\lim_{\kk \to \infty} \omega\args{ \vxNk; \vf, \linFeasApproxk, \normTwo{\bullet}} = 0.
\]
\end{lemm}

\begin{proof}
Similar to~\cite[Lemma~7]{berkemeierDerivativeFreeMultiobjectiveTrust2021}, we can show that for $\kk\in \mathbb{N}_0$ it must hold that
\[
\abs{
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normTwo{\bullet}} -
\omega\args{ \vxNk; \vf, \linFeasApproxk, \normTwo{\bullet}}} =
\abs{
\omegak_2 - \omegak_2}
\le \mathtt{s} \Deltak
\]
for some constant $\mathtt{s} > 0$ (independent of $\kk$). The triangle inequality yields
\[
\abs{\omegak_2}
\le
\abs{\omegak_2 - \omegak_2} + \abs{\omegak_2} \le
\mathtt{s}\Deltak + \abs{\omegak_2}
\]
and the \ac{rhs} goes to zero for $k\to \infty$.
\end{proof}

Note that \cref{thm:true_objective_criticality} applies both to a sequence of the main algorithm and to an infinite subsequence of the Criticality Routine (see \cref{thm:infinite_crit_routine}). In our second step, we now replace the \emph{constraint} surrogates with the original functions.

\begin{lemm}\label{thm:extinesHeroizes}
Suppose \cref{ass:covering,ass:lipschitz_gradients_true,ass:fully_linear_models,ass:normal_step_exists_bounded} hold. Let $\{\vxk\}\subseteq \mathbb{R}^\varDim$ be an algorithmic sequence with $\vxk\to \bar{\vx} \in \feas$, $\omega\args{ \vxNk; \vf, \linFeasApproxk, \normInf{\bullet}} \to 0$ and $\Deltak \to 0$. Further, assume that the \acp{mfcq} hold at $\bar{\vx}$. i.e., the rows of $\DH(\bar{\vx})$ are linearly independent, and there is a direction $\vd \in \mathbb{R}^\varDim$ such that $\DH(\bar{\vx})\vd = \vZ$ and $\noDotProd{\vd}{\dg_\indg(\bar{\vx})} < 0$ for all $\indg$ with $\gL_\indg(\bar{\vx}) = 0$. Then $\bar{\vx}$ is a KKT-point of~\eqref{eqn:mop_nonlin}.
\end{lemm}

\begin{proof}
In the following, we use $\DmH$ and $\DmG$ to denote the Jacobians of the surrogates $\mhk$ and $\mgk$, evaluated at $\vxk$. Likewise, $\DFk$ is the Jacobian of the true objective function at $\vxk$, while $\mgkSub \le \vZ_\gDim$ is defined as $\mgk(\vxk) + \DmG \cdot \vnk$.

\Cref{ass:normal_step_exists_bounded} gives $\vnk \to \vZ$. Because of this, and the error bounds in \cref{ass:fully_linear_models}, we have that $\vxk \to \bar{\vx}$ as $\Deltak \to 0$, as well as $\DmH \to \DH = \jac \vh (\bar{\vx})$, $\DmG\to \DG = \jac \vh (\bar{\vx})$ and $\mgkSub\to \vg = \vg(\bar{\vx})$.

By assumption, $\omega\args{ \vxNk; \vf, \linFeasApproxk, \normInf{\bullet}} \to 0$. We thus have a sequence of linear programs,
\begin{equation*}
\omega\args{ \vxNk; \vf, \linFeasApproxk, \normInf{\bullet}} =
\max_{\vd \in \mathbb R^n, \beta^- \in \mathbb R}
\begin{bmatrix}
\tran{\vZ}_{\varDim} & 1
\end{bmatrix}
\begin{bmatrix}
\vd \\ \beta^-
\end{bmatrix}
\quad\text{s.t.}\quad
\begin{bmatrix}
-\eye{\varDim} & \vZ_\varDim \\
\eye{\varDim} & \vZ_\varDim \\
\DFk & \ve{1}_\fDim \\
\DmH & \vZ_\hDim \\
\DmG & \vZ_\gDim
\end{bmatrix}
\begin{bmatrix}
\vd \\ \beta^-
\end{bmatrix}
\begin{matrix}
\le \\ \le \\ \le \\ = \\ \le
\end{matrix}
\begin{bmatrix}
\ve{1}_\varDim\\
\ve{1}_\varDim\\
\vZ_\fDim \\
\vZ_\hDim\\
-\mgkSub \\
\end{bmatrix},
\end{equation*}
the values of which go to zero. The dual problems are
\begin{equation}
\begin{aligned}
&\min_{
\ve{y}^1, \ve{y}^2, \ve{y}^3, \ve{y}^5 \ge \vZ, \ve{y}^4 \in \mathbb R^\hDim}
\begin{bmatrix}
\tran{\ve{1}}_{\varDim} & \tran{\ve{1}}_{\varDim} & \tran{\vZ}_{\fDim} & \tran{\vZ}_{\hDim} & - \tran{\mgkSub}
\end{bmatrix}
\begin{bmatrix}
\ve{y}^1 \\ \vdots \\ \ve{y}^5
\end{bmatrix}
\quad \text{s.t.}\quad
\begin{bmatrix}
-\eye{\varDim} & \eye{\varDim} & \tran{\DFk} & \tran{\DmH} & \tran{\DmG}\\
\tran{\vZ}_{\varDim} & \tran{\vZ}_{\varDim} & \tran{\ve{1}}_{\fDim} & \tran{\vZ}_{\hDim} & \tran{\vZ}_{\gDim}
\end{bmatrix}
\begin{bmatrix}
\ve{y}^1 \\ \vdots \\ \ve{y}^5
\end{bmatrix} =
\begin{bmatrix}
\vZ_{\varDim} \\ 1
\end{bmatrix}.
\end{aligned}\label{eqn:dual_k_unscaled}
\tag{D\ensuremath{_k}}
\end{equation}
By strong duality, the dual problem for $k\in \mathbb{N}_0$ is also always feasible, and its optimal value equals the primal optimal value. Suppose that $\set{ (\ve{y}_k^1, \ldots, \ve{y}_k^5)}$ is a sequence of dual optimizers. By strong duality, it follows from $\omega\args{ \vxNk; \vf, \linFeasApproxk, \normInf{\bullet}} \to 0$ that
\begin{equation*}
\lim_{k\to \infty} \ve{y}_k^1 = \vZ,
\quad
\lim_{k\to \infty} \ve{y}_k^2 = \vZ,
\quad\text{and}\quad
\lim_{k\to \infty} - \noDotProd{\mgkSub}{\ve{y}_k^5} = 0.
\end{equation*}
First, we show that it must also hold for $i\in \{3,4,5\}$ that the sequences $\{ \ve{y}^i\}$ are bounded (see also~\cite{echebestInexactRestorationDerivativefree2017} for a similar idea). To derive a contradiction, assume that for some $i \in \{ 3, 4, 5\}$ the sequence $\{ \ve{y}_k^i\}$ is not bounded. Define
\(
\nu_k = \max \set{ 1, \max_{i=3,4,5} \normInf{\ve{y}_k^i}},
\) and we get a bounded sequence of scaled variables
\[
\tilde{y}_k^i = \frac{1}{\nu_k} y_k^i \qquad i = 1,\ldots,5.
\]
Then, $\nu_k \xrightarrow{k\to \infty} \infty$, but we can take a subsequence $\mathcal K$ of indices such that for every $i\in\{3,4,5\}$ the sequence ${\{\tilde{\ve{y}}_k^i\}}_{k\in \mathcal K}$ converges to $\bar{\ve{y}}^i$ and so that for one $i \in \{3,4,5\}$ the limit $\bar{\ve{y}}^i$ is not zero. Consequently, we have
\begin{equation}
\lim_{\substack{k\to\infty,\\ k\in \mathcal K}} (\tilde{\ve{y}}_k^1,\ldots, \tilde{\ve{y}}_k^5) = (\bar{\ve{y}}^1, \ldots, \bar{\ve{y}}^5) = (\vZ_{\varDim}, \vZ_{\varDim}, \bar{\ve{y}}^3, \bar{\ve{y}}^4, \bar{\ve{y}}^5),
\quad\text{and}\quad
\lim_{\substack{k\to\infty,\\ k\in \mathcal K}} - \noDotProd{\mgkSub}{\tilde{\ve{y}}_k^5} = -\noDotProd{\ve{g}}{\ve{\bar{y}}^5} = 0.\label{eqn:week_afford_income}
\end{equation}
For any dual feasible $[\ve{y}_k^1, \ldots, \ve{y}_k^5]$, the vector $[\tilde{\ve{y}}_k^1, \ldots, \tilde{\ve{y}}_k^5]$ is feasible for the problem
\begin{equation*}
\begin{aligned}[t]
&
\min_{
\tilde{\ve{y}}^1,
\tilde{\ve{y}}^2,
\tilde{\ve{y}}^3,
\tilde{\ve{y}}^5 \ge \vZ,
\tilde{\ve{y}}^4 \in \mathbb R^\hDim}
\begin{bmatrix}
\tran{\ve{1}}_{\varDim} &
\tran{\ve{1}}_{\varDim} &
\tran{\vZ}_{\fDim} &
\tran{\vZ}_{\hDim} & -
\tran{\mgkSub}
\end{bmatrix}
\begin{bmatrix}
\tilde{\ve{y}}^1 \\ \vdots \\ \tilde{\ve{y}}^5
\end{bmatrix}
\quad \text{s.t.}\quad
\begin{bmatrix}
\eye{\varDim} &
\eye{\varDim} & \tran{\DFk} & \tran{\DmH} & \tran{\DmG} \\
\tran{\vZ}_{\varDim} & \tran{\vZ}_{\varDim} & \tran{\ve{1}}_{\fDim} & \tran{\vZ}_{\hDim} & \tran{\vZ}_{\gDim}
\end{bmatrix}
\begin{bmatrix}
\tilde{\ve{y}}^1 \\ \vdots \\ \tilde{\ve{y}}^5
\end{bmatrix} =
\begin{bmatrix}
\vZ_{\varDim} \\ \nicefrac{1}{{\nu}_k}
\end{bmatrix}.
\end{aligned}
\end{equation*}
In the limit, the last constraint becomes $\dotProd{{\ve{1}}_{\fDim}}{\bar{\ve{y}}^3} = 0$, the limiting problem is
\begin{equation}
\begin{aligned}[t]
&
\min_{
\tilde{\ve{y}}^1,
\tilde{\ve{y}}^2,
\tilde{\ve{y}}^3,
\tilde{\ve{y}}^5 \ge \vZ,
\tilde{\ve{y}}^4 \in \mathbb R^\hDim}
\begin{bmatrix}
\tran{\ve{1}}_{\varDim} &
\tran{\ve{1}}_{\varDim} &
\tran{\vZ}_{\fDim} &
\tran{\vZ}_{\hDim} & -\tran{\vg}
\end{bmatrix}
\begin{bmatrix}
\tilde{\ve{y}}^1 \\ \vdots \\ \tilde{\ve{y}}^5
\end{bmatrix}
\quad \text{s.t.}\quad
\begin{bmatrix}
\eye{\varDim} &
\eye{\varDim} & \tran{\DF} & \tran{\DH} & \tran{\DG} \\
\tran{\vZ}_{\varDim} & \tran{\vZ}_{\varDim} & \tran{\ve{1}}_{\fDim} & \tran{\vZ}_{\hDim} & \tran{\vZ}_{\gDim}
\end{bmatrix}
\begin{bmatrix}
\tilde{\ve{y}}^1 \\ \vdots \\ \tilde{\ve{y}}^5
\end{bmatrix} =
\begin{bmatrix}
\vZ_{\varDim} \\ 0
\end{bmatrix}.
\end{aligned}
\tag{D$'$}\label{eqn:dual_scaled}
\end{equation}
The primal thus has constant objective value:
\begin{equation}
\begin{aligned}[t]
\max_{\vd \in \mathbb{R}^\varDim, \beta^- \in \mathbb R} 0
\qquad\text{s.t.}\qquad
\begin{bmatrix}
-\eye{\varDim} & \vZ_{\varDim} \\
\eye{\varDim} & \vZ_{\varDim} \\
\DF & \ve{1}_{\fDim} \\
\DH & \vZ_{\hDim} \\
\DG & \vZ_{\gDim}
\end{bmatrix}
\begin{bmatrix}
{\vd} \\ {{\beta}}^-
\end{bmatrix}
\begin{matrix}
\le \\ \le \\ \le \\ = \\ \le
\end{matrix}
\begin{bmatrix}
\ve{1}_{\varDim} \\ \ve{1}_{\varDim} \\ \vZ_{\fDim} \\ \vZ_{\hDim} \\ - \ve{g}
\end{bmatrix}.
\end{aligned}\label{eqn:primal_scaled}
\tag{P$'$}
\end{equation}
By upper semi-continuity of the feasible set mappings (see~\cite{wetsContinuityValueLinear1985a}), $\bar{\ve{y}}$ is feasible for~\eqref{eqn:dual_scaled}, and we see from~\eqref{eqn:week_afford_income} and~\eqref{eqn:primal_scaled} that it is also optimal. Let $\vd\in \mathbb R^n$ be a direction with $\normInf{\vd} \leq 1$ adhering to the
\acp{mfcq}. Furthermore, let $\beta$ be such that $(\vd, \beta)$ is feasible for~\eqref{eqn:primal_scaled}. Then $(\vd, \beta)$ is optimal for~\eqref{eqn:primal_scaled}, and we make the following observations:
\begin{itemize}
\item \textbf{$\bar{\ve{y}}^3$ must be zero,} because of the second constraint in~\eqref{eqn:dual_scaled}.
\item \textbf{$\bar{\ve{y}}^5$ must be zero:} Because of the \acp{mfcq} it holds that $\DG_{[\indg, :]} \cdot \vd < 0$ whenever $\gL_\indg = 0$. By complementary slackness, it follows for these indices $\indg$ that the entries in $\bar{\ve{y}}^5$ must be $0$. From~\eqref{eqn:week_afford_income} it also follows for the other indices that their values must be $0$.
\item But with $\bar{\ve{y}}^4 \ne \vZ$ it then follows from the first constraint in~\eqref{eqn:dual_scaled} that $\dotProd{\DH}{\bar{\ve{y}}^4} = \vZ$ in contradiction to the \acp{mfcq}.
\end{itemize}
Hence, for all $i \in \{3,4,5\}$ the sequence $\{\ve{y}_k^i\}$ of (unscaled) Lagrange multipliers must also be bounded!

We thus can take a subsequence $\mathcal K$ so that ${\{\ve{y}_k\}}_{k\in \mathcal K}$ converges to some $(\vZ, \vZ, \bar {\ve{y}}_k^3, \bar{\ve{y}}_k^4, \bar{\ve{y}}^5)$ and (by upper semi-continuity of the feasible set mapping) the limit point is feasible for the limiting problem of~\eqref{eqn:dual_k_unscaled}, which happens to be~\eqref{eqn:buntingGeologically}. The optimal value at $(\vZ, \vZ, \bar {\ve{y}}_k^3, \bar{\ve{y}}_k^4, \bar{\ve{y}}_k^5)$ is 0 due to strong duality. Because the corresponding primal is~\eqref{eqn:legitimizersSannop}, it follows from \cref{thm:subproblems_implies_kkt} that $\bar{\vx}$ is a KKT point of~\eqref{eqn:mop_nonlin}.
\end{proof}

With~\cref{thm:subproblems_implies_kkt}, it is now easy to derive the main result:
\begin{theo}[Convergence to KKT-points]
Suppose the same assumptions as in~\cref{thm:convergence1} hold and that $\{\vxk\}$ is a quasi-stationary subsequence with limit point $\bar{\vx}$. If the \acp{mfcq} hold at $\bar{\vx}$, then $\bar{\vx}$ is a KKT-point of~\eqref{eqn:mop_nonlin}.
\end{theo}

\section{Numerical Examples}\label{sec:examples}

In this section, we provide numerical examples in order to evaluate the strengths and weaknesses of our algorithm. There are two examples using an implementation close to the description in~\cref{section:algorithm}. Additionally, we give comparisons with the \dfmo{} algorithm~\cite{liuzziDerivativeFreeApproachConstrained2016} on the \testmo{} set of test problems. This comparison is done not with the main algorithm, but with algorithms derived from it in order to generate multiple solutions.

\subsection{Implementation Details}

All algorithms are available in the research package \compromise{}~\cite{berkemeierCoMPrOMISE2024} in the Julia language. The code to generate the experiments can be found in a separate GitHub repository.\footnote{\url{https://github.com/manuelbb-upb/MOBenchmarks}} It also contains utilities to interact with the Fortran implementations of \dfmo{} and
\testmo{} provided by the authors of~\cite{liuzziDerivativeFreeApproachConstrained2016}.

Trust-regions are boxes, but we allow for the use of other norms in the objective of~\eqref{eqn:itr_normal} and the constraint in~\eqref{eqn:itr_tangential}. In fact, by default linear norms are used because this appears to improve the quality of solutions with the inner solver \texttt{HiGHS}~\cite{huangfuParallelizingDualRevised2018}. For restoration, we currently use \texttt{COBYLA}~\cite{powellDirectSearchOptimization1994} from the \texttt{NLopt} package~\cite{johnsonNLoptNonlinearoptimizationPackage2007}.

The version of \compromise{} used in this demonstration allows for \ac{rbf} surrogate models according to the construction in~\cite{wildDerivativeFreeOptimizationAlgorithms2009}. It is also possible to use first or second degree Taylor models. These need derivative information, which we provide by means of finite difference approximations using \texttt{FiniteDiff.jl}.

The plots are generated with \texttt{Makie.jl}~\cite{danischMakiejlFlexibleHighperformance2021}.

\subsection{Constrained Two-Parabolas Problem}

Variations of the two parabolas problem are popular test cases for \ac{moo} methods. The objectives are simply two $\varDim$-variate quadratic functions and (in the unconstrained case) the Pareto Set is the line connecting their respective minima.
\subsubsection{Inequality Constraint} The following constrained version is taken from~\cite{gebkenDescentMethodEquality2018}.
\begin{equation}
\begin{aligned}
&\min_{\vx\in \mathbb{R}^2}
\begin{bmatrix}
(x_1 - 2)^2 + (x_2 - 1)^2 \\
(x_1 - 2)^2 + (x_2 + 1)^2
\end{bmatrix}
\quad\text{s.t.}\quad
\gL(\vx) = 1 - x_1^2 - x_2^2 \le 0.
\end{aligned}\label{eqn:bennets_problem}
\tag{TPineq}
\end{equation}
The feasible set of~\eqref{eqn:bennets_problem} is $\mathbb{R}^2$ without the interior of the $\normTwo{\bullet}$ unit ball, so it is non-convex. The Pareto critical set is the line connecting $[2, -1]$ and $[2, 1]$ and the left boundary of the unit ball:
\begin{equation*}
\mathcal P_c =
\set{
\begin{bmatrix}
2\\
s
\end{bmatrix}: s\in[-1,1]}
\bigcup
\set{
\begin{bmatrix}
\cos t\\
\sin t
\end{bmatrix} : t \in [\pi - \theta, \pi + \theta], \theta = \arctan\args{\frac{1}{2}}}.
\end{equation*}
We have plotted the critical set with black lines in~\cref{fig:bennets_problem_ineq}.

\begin{figure}[h]
\centering
\includegraphics[width=\textwidth]{imgs/two_paraboloids_ineq.pdf}
\caption{Optimization trajectories with different surrogate types for~\eqref{eqn:bennets_problem}. The upper legend shows the number of function calls in parentheses. Different strokecolors correspond to different values of the threshold parameter $\CCsuccess$.}\label{fig:bennets_problem_ineq}
\end{figure}

For this problem, we first apply our algorithm three times with the same parameters and beginning at $\vx_0 = [-2, 0.5]$, but with different model types. We compare \ac{rbf} surrogate models with first and second degree Taylor polynomials of the objective and constraint functions. The \ac{rbf} models use the cubic radial function $r \mapsto r^3$. The derivatives for the Taylor models are approximated using finite differences (to conform to the assumption that gradients are not available exactly). The algorithm parameters are
\begin{equation}
\begin{gathered}
\DeltaInit =.5,
\DeltaMax = 16,
\gammass = 0.25,
\gammas = 0.5,
\gammag = 2.0,
\filterConst = 0.1,
\CCaccept = 0.01,
\CCcritTestchi = 10^{-4},
\CCcritTesttheta = 10^{-6},\\
\kappa_\theta = 10^{-4},
\psi = 2,
\CCBeta = 1000,
\CCMu = 3000,
\CCalpha = 0.5,
\CCDelta = 0.99,
\CCmu = 100,
\mu = 0.01.
\end{gathered}\label{eqn:kafilaDisgowned}
\end{equation}
The parameter $\CCdeltan$ is not required for computations, but assumed to have a value $\CCdeltan \ge \CCcritTesttheta$. These settings largely follow the recommendations in~\cite{fletcherGlobalConvergenceTrustRegion2002}. For these runs, the success threshold is $\CCsuccess = 0.45$. We stop after 100 iterations or if the trust region radius is reduced to below $10^{-9}$. Afterwards, we change the starting point to $x_0 = [-2, -0.5]$ and the success threshold to $\CCsuccess = 0.9$ for three runs with the same kinds of surrogates. Finally, there is also a run starting at $x_0 = [-2, 0]$ with $\CCsuccess =.45$.

The results are depicted in~\cref{fig:bennets_problem_ineq}. As can be seen, all runs converge to a point on the Pareto critical set and avoid the infeasible area. In terms of objective evaluations, the \ac{rbf} models need the fewest function calls. In terms of iterations, the second degree Taylor polynomials perform best.

Almost all solutions lie on the connecting line of the individual objective minima. The only run converging to the left part of the critical set started out close to it. Of course, this behavior is in part caused by the problem geometry and the surrogate construction methods. But if our wish is to find non-dominated points rather than the superset of critical points, we benefit from the descent nature of the method.

Compared to the Taylor models, the RBF models sometimes take large steps that do not necessarily reduce both objectives to the same degree. This is likely due to them providing more inexact gradients. Moreover, using $\normInf{\bullet}$ in~\eqref{eqn:itr_tangential} induces a bias towards the corners of the unit ball.

We can also see the zig-zag motion close to critical points that is typical for gradient-based algorithm. It is very pronounced for second degree Taylor models. The Taylor models of degree~2 can be considered almost exact for~\eqref{eqn:bennets_problem}. Hence, the trial points achieve the predicted objective decrease and nearly every iteration is deemed successful, which makes the trust-regions have large radius and practically reduces our method to ``vanilla'' gradient-descent.

Finally, it is worthwile to compare the iteration trajectories to those in~\cite{gebkenDescentMethodEquality2018}, where an active set strategy is used. In contrast to the results presented there, we stay further away from the infeasible area. A possible explanation lies in the fact that (in our algorithm) steps are computed with inexact approximations of \emph{every} constraint function, instead of only considering the active constraints.

\subsubsection{Equality Constraint} We can take the same objectives and turn the constraint into an equality constraint:
\begin{equation}
\begin{aligned}
&\min_{\vx\in \mathbb{R}^2}
\begin{bmatrix}
(x_1 - 2)^2 + (x_2 - 1)^2 \\
(x_1 - 2)^2 + (x_2 + 1)^2
\end{bmatrix}
\quad\text{s.t.}\quad
\gL(\vx) = 1 - x_1^2 - x_2^2 = 0.
\end{aligned}\label{eqn:bennets_problem_eq}
\tag{TPeq}
\end{equation}
Now, the feasible set is the $\normTwo{\bullet}$ unit sphere. Subsequently, the connecting line of individual objective minima cannot be optimal anymore. Analysis shows that the Pareto critical set now consists of two segments of the unit sphere:
\begin{equation*}
\mathcal P_c =
\set{
\begin{bmatrix}
\cos t\\
\sin t
\end{bmatrix} : t \in [\pi - \theta, \pi + \theta]
\cup [2\pi - \theta, 2\pi + \theta],\;
\theta = \arctan\args{\frac{1}{2}}}.
\end{equation*}
We mostly keep the parameters~\eqref{eqn:kafilaDisgowned} and set $\CCsuccess = 0.9$. This time, for different starting points, we vary the initial trust region radius and the compatibility parameter $\CCDelta$. Doing so allows us to illustrate the effect of these parameters on the restoration procedure (amongst others) by comparing trajectories and pairs of trajectories.

The results are depicted in~\cref{fig:bennets_problem_eq}. First, this example illustrates well how restoration steps return the iterations back close to the feasible set. Secondly, we again notice a tendency towards the objective minima; only the green trajectory ends on the left part of the Pareto Set. This is due to the significantly smaller initial trust region radius. Comparing the blue trajectory with the others confirms that a greater $\CCDelta$ value is more permissive with respect to declaring the normal step compatible. Specifically, the purple trajectory has been set up symmetrically to the blue one. But due to $\CCDelta$ being smaller, purple starts with a restoration right away. Orange starts near blue and can avoid initial restoration because it starts closer to the feasible set. But as with purple there are noticeably more restoration iterations overall, compared to blue. Comparing orange and green highlights that \enquote{compatibility} is a function not only of $\CCDelta$ but of $\Deltak$ as well. That is why for green, with the smaller initial radius, the first iteration is a restoration.

\begin{figure}[h]
\centering
\includegraphics{imgs/two_paraboloids_eq.pdf}
\caption{Four optimization trajectories for~\eqref{eqn:bennets_problem_eq} with cubic RBF surrogates and varying initial trust region radii and compatibility parameters.}\label{fig:bennets_problem_eq}
\end{figure}

\subsection{Non-Convex Test-Problem ``W3''}\label{sec:experiments_w3}

\begin{figure}[h]
\centering
\includegraphics{imgs/mw3_trajectories.pdf}
\caption{Results for the MW3 problem and comparison with the weighted-sum approach. Different marker symbols refer to different optimization techniques. Each initial objective vector is connected to the solution vector by a thin line. For \compromise{}, the lines represent the image space trajectories. The gray area show all objective vectors that are feasible if we read the constraints as functions of $f_1$ and $f_2$. The darker area actually has valid pre-images.}\label{fig:mw3}
\end{figure}

The second example is problem ``W3'' taken from~\cite{maEvolutionaryConstrainedMultiobjective2019}. The problem has box constraints and two non-linear inequality constraints. Whilst the unconstrained Pareto Front is simply a line in objective space, the constraints make it partially non-convex. We have chosen the problem with $3$ variables and $2$ objectives to compare our algorithm against a simple weighted-sum scalarization. It reads:
\begin{equation}
\begin{aligned}
&\min_{\vx \in [0, 1]^3}
\begin{bmatrix}
f_1(\vx) \\
f_2(\vx)
\end{bmatrix} =
\min_{\vx \in [0, 1]^3}
\begin{bmatrix}
x_1\\
d(\vx)\cdot\args{1 - \frac{f_1(\vx)}{d(\vx)}}
\end{bmatrix}
\qquad&\text{s.t.}
\\
&c_1(\vx) = f_1(\vx) + f_2(\vx) - 1.05 - 0.45\sin(0.75 \pi \cdot l(\vx))^6 \le 0, \\
&c_2(\vx) = -f_1(\vx) - f_2(\vx) + 0.85 + 0.3\sin(0.75 \pi \cdot l(\vx))^2 \le 0,\\
&d(\vx) = 1 + 2(x_2 + (x_1 - 0.5)^2 - 1)^2 +2(x_3 + (x_2 - 0.5)^2 - 1)^2, \\
&l(\vx) = \sqrt{2}(f_2(\vx) - f_1(\vx)).\\
\end{aligned}
\tag{MW3}\label{eqn:mw3}
\end{equation}
Because the constraints $c_1$ and $c_2$ describe the feasible set in terms of the objective values we can show the attainable objective values in objective space in~\cref{fig:mw3}. The critical set is again shown in black.

For 15 fixed starting points, we apply 3 different algorithms:
\begin{itemize}
\item Our trust-region algorithm \compromise{} with default settings,
\item a modified version of \compromise{}, with normalized gradients and requiring strict descent,
\item and a weighted-sum minimization of $\nicefrac{1}{2} f_1 + \nicefrac{1}{2} f_2$ with \texttt{COBYLA}.
\end{itemize}

As we can see from the results in~\cref{fig:mw3}, the weighted-sum approach is only able to find convex parts of the Pareto Front. Both multi-objective descent techniques are able to find solutions on the non-convex parts. The standard version of \compromise{} has only one solution in the non-convex region. In the standard version, we allow partial descent. By requiring strict descent, the modified version does not “run away” from the non-convex parts. Additionally, normalizing gradients before computing the steepest descent direction centers the steepest descent direction. Nonetheless, both versions of \compromise{} show a typical phenomenon of gradient-based multi-objective optimization techniques: Some solutions are on the boundary and only weakly optimal.

\subsection{Comparison with \dfmo{} on \testmo{} Suite}

\subsubsection{Solvers}

The solver \dfmo{} is a Fortran implementation of a derivative-free multi-objective optimization algorithm for problems with box constraints and nonlinear inequality constraints~\cite{liuzziDerivativeFreeApproachConstrained2016}. The code is available under an open license\footnote{\url{https://github.com/DerivativeFreeLibrary/DFMO}}, which makes it a convenient candidate to test our algorithm against.
\dfmo{} combines a penalty approach with inexact line-search along quasi-random and coordinate directions to find a set of critical solutions. However, our method produces only a single critical output. We thus designed various \enquote{outer} algorithms based on the same principles.

The \emph{sequential} outer algorithm starts with a set of initial points and sequentially performs inner updates according to the main algorithm from \cref{section:algorithm}. That is, for each point, first a surrogate is constructed, and then we compute an inexact normal step and an inexact descent step. It is very similar to just executing multiple runs with the main algorithm, but allows us to share the evaluation database in a more effective manner. We refer to the sequential outer algorithm as \cseq{}. A major drawback of \cseq{} becomes obvious if nonlinear constraints are present: A single restoration phase can easily use up any computational budget.

That is why we also test \cpen{}, the sequential outer algorithm applied to the box-constrained problem with penalty-augmented objectives
\[
Z_\indf(\vx) = \fL_\indf(\vx) +
\frac{1}{\varepsilon}
\sum_{\indg=1}^\gDim \max \set{0, \gL_\indg(\vx)}
\qquad\text{for some $\varepsilon > 0$.}
\]
The experiments are performed with $\varepsilon = 0.1$. Note that \cpen{} does no longer satisfy the assumption that the objectives are smooth. In between outer iterations we may perform non-dominance testing to not waste time on dominated solutions.

With \cpen{} we abandon the composite step approach, as no normal steps are required. But it seems reasonable that normal steps are beneficial to finding feasible solutions by exploiting the local structure of the feasible set. For this reason, we thought of a \emph{set-based} outer algorithm, \cset{}. In contrast to \cseq{}, more information between solutions is shared and we try to avoid restoration as much as possible. A high-level description of \cset{} is given in the Appendix, \cref{sec:outerAlgo}. Unfortunately, there is no convergence theory for it yet.

Compared to \dfmo{}, all our algorithms suffer from a lack of \emph{exploration} mechanisms. Any exploration happens by accident due to inexact surrogates. We try to compensate by starting with large initial populations.

\subsubsection{Setup and Metrics}\label{sec:setup_and_metrics}
We use the test problem library \testmo{} and largely follow the setup in~\cite{liuzziDerivativeFreeApproachConstrained2016}. From \testmo{}, we take those 51 multi-objective problems with 3 or more variables. All problems have finite box constraints. We also generate 214 test problems with non-linear constraints by augmenting each of the 51 problems with one of six constraint functions provided in~\cite{liuzziDerivativeFreeApproachConstrained2016} and discarding infeasible problems. Each solver has a budget of 2000 evaluations.
\dfmo{} is run with default settings. The algorithms based on \compromise{} have small acceptance threshold $\CCaccept = 10^{-4}$ but a large success threshold $\CCsuccess = 0.9$. All functions are modelled with cubic \ac{rbf} surrogates. For the \compromise{} algorithms, we use the first $\min \set{500, 100\varDim}$ Halton points as starting points.

We present the results of our experiments by means of performance profiles~\cite{custodioDirectMultisearchMultiobjective2011, dolanBenchmarkingOptimizationSoftware2002} for performance indicators derived from different metrics. Loosely speaking, we define a performance indicator as a value $t_{p, s}$ that is lower the better solver $s \in \solvers$ performs on problem $p \in \problems$. The values are then transformed to obtain normalized performance values $r_{p,s}$. Finally, we plot the distribution functions
\[
\rho_s(\tau) = \frac{1}{\abs{\problems}} \abs{\set{p \in \problems: r_{p,s} \leq \tau}},
\]
starting at the lowest possible $\tau$. Usually, $r_{p,s}$ is the performance \emph{ratio}, obtained by dividing $t_{p,s}$ by the best possible value obtained among all solvers, i.e.,
\begin{equation}
r_{p,s} = \frac{t_{p,s}}{\min\set{t_{p, s}:\; s \in \solvers}}.
\tag{T1}\label{eq:ratio_transform}
\end{equation}
Then, $\tau \geq 1$ and $\rho_s(1)$ is the percentage of problem instances for which $s$ was the best. The solver $s$ with highest $\rho_s(1)$ value is the best overall, whilst high values $\rho_s(\tau)$ for $\tau > 1$ indicate solver robustness. However, if $t_{p,s}$ can become $0$, we instead use a min-max-normalization, i.e.,
\begin{equation}
r_{p,s} =
\begin{cases}
0 & \text{if $\bar{t}_p - \underline{t}_p = 0$},\\
\frac{t_{p,s} - \underline{t}_p}{\bar{t}_p - \underline{t}_p} & \text{otherwise},
\end{cases}
\quad
\underline{t}_p = \min\set{t_{p, s}:\; s \in \solvers},
\bar{t}_p = \max\set{t_{p, s}:\; s \in \solvers}.
\tag{T2}\label{eq:minmax_transform}
\end{equation}
Subsequently, $r_{p,s} \in [0, 1]$ and $\tau\in [0, 1]$, and $r_{p,s}$ relates the performance of $s$ to both the best and worst solver for problem instance $p$.

From~\cite{custodioDirectMultisearchMultiobjective2011} we take the purity metric, and the spread metrics $\spreadGamma$ and $\spreadDelta$. All these metrics require a reference front $\front_p$ or extreme points. We build the reference front $\front_p$ by taking the union of the approximate fronts $\front_{p,s}$ returned by all solvers and removing the dominated points from this set.

The purity of a solver on a problem is the percentage of solutions of this solver that are also in the reference front. To make it suitable for performance profiles, $t_{p,s}$ is the reciprocal of that fraction, i.e., $t_{p,s} = \nicefrac{\abs{\front_p}}{\abs{\front_{p,s}\cap \front_p}}$, and $r_{p,s}$ is obtained from~\eqref{eq:ratio_transform}.

The $\spreadGamma$ metric relates to the maximum inter-solution distance within front approximations. The metric $\spreadDelta$ is an indicator for how well approximate solution value vectors are distributed. In both cases $t_{p,s}$ can become zero, thus~\eqref{eq:minmax_transform} is used.

We consider two more metrics. The hyperarea ratio~\cite{audetPerformanceIndicatorsMultiobjective2020} is an indicator of how much volume is dominated by a solver front compared to the reference front. Like with the purity, we take the reciprocal and use~\eqref{eq:ratio_transform}. The reference point for the hypervolume computation is constructed like in~\cite{brasUsePolynomialModels2020}.

Moreover, we take the median of the constraint violation values of approximate solutions as a performance indicator, as well as the 3rd quartile values. In this case, we measure constraint violation as $\sum_\indg\max \set{0, \gL_\indg(\vx)}$. These values are normalized using~\eqref{eq:minmax_transform}.

For every indicator except \enquote{constraint violation}, we only take into consideration those solutions where the constraint violation is below $10^{-3}$.

\subsubsection{Results and Discussion}

\begin{figure}[h]
\centering
\includegraphics[width=0.32\textwidth]{imgs/purity_DFMO_CSET_unconstr.pdf}
\includegraphics[width=0.32\textwidth]{imgs/purity_DFMO_CSEQ_unconstr.pdf}
\includegraphics[width=0.32\textwidth]{imgs/purity_CSEQ_CSET_unconstr.pdf}
\caption{Pairwise purity comparison on unconstrained problems with $\tau \ge 1$.}\label{fig:purity_unconstr}
\end{figure}

\begin{figure}[h]
\centering
\includegraphics[width=0.32\textwidth]{imgs/purity_DFMO_CSET_constr.pdf}
\includegraphics[width=0.32\textwidth]{imgs/purity_DFMO_CSEQ_constr.pdf}
\includegraphics[width=0.32\textwidth]{imgs/purity_DFMO_CPEN_constr.pdf}
\\%
\includegraphics[width=0.32\textwidth]{imgs/purity_CSEQ_CSET_constr.pdf}
\includegraphics[width=0.32\textwidth]{imgs/purity_CSEQ_CPEN_constr.pdf}
\includegraphics[width=0.32\textwidth]{imgs/purity_CSET_CPEN_constr.pdf}
\caption{Pairwise purity comparison on constrained problems with $\tau \ge 1$.}\label{fig:purity_constr}
\end{figure}

As can be seen in \cref{fig:purity_unconstr} and \cref{fig:purity_constr}, \dfmo{} consistently outperforms all other solvers. The sequential outer algorithms seem to do very bad. Fortunately, \cset{} appears to find solutions that have better quality. It is still outperformed by \dfmo{}, which might be due to \cset{} being more prone to local minima.

The plots for the hyperarea ratio indicator in \cref{fig:hv} confirm that \dfmo{} produces the best solutions with \cset{} as the runner-up.
\begin{figure}[h]
\centering
\includegraphics[width=0.48\textwidth]{imgs/inv_hyperarea_ratio_unconstr.pdf}
\includegraphics[width=0.48\textwidth]{imgs/inv_hyperarea_ratio_constr.pdf}
\caption{Performance based on the inverse hyperarea ratio for constrained (c) and unconstrained (u) problems with $\tau \ge 1$.}\label{fig:hv}
\end{figure}

A similar trend is seen with the $\spreadGamma$ spread in \cref{fig:spread}, both in the constrained and unconstrained case. Again, \dfmo{} outperforms all other solvers. Due to the quasi-random directions it uses, and its expansive line-search routine, it has much better exploration capabilities than the \compromise{} algorithms.
\\
At first glance, the sequential outer algorithms \cseq{} and \cpen{} seem to do better regarding the $\spreadDelta$ spread. But this is a fluke. The sequential algorithms usually find way fewer (and worse) solutions, whilst \dfmo{} and \cset{} obtain more solutions that are sometimes clustered, with gaps in-between clusters. These gaps result in high $\spreadDelta$ values, even if the cluster-gaps for \dfmo{} or \cset{} are smaller than inter-solution gaps for \cseq{} or \cpen{}.
\begin{figure}[h]
\centering
\includegraphics[width = 0.48\textwidth]{imgs/xi_spread_unconstr.pdf}
\includegraphics[width = 0.48\textwidth]{imgs/xi_spread_constr.pdf}
\\%
\includegraphics[width = 0.48\textwidth]{imgs/theta_spread_unconstr.pdf}
\includegraphics[width = 0.48\textwidth]{imgs/theta_spread_constr.pdf}
\caption{Spread performance for constrained (c) and unconstrained (u) problems with $\tau \in [0, 1]$.}\label{fig:spread}
\end{figure}

Finally, both \dfmo{} and \cset{} find mostly feasible solutions, see \cref{fig:constraint_viol}. Surprisingly, \cpen{} is not too bad either. On the other hand, \cseq{} behaves as expected. Likely, much of the budget is spent on restoration a single solution, or very few solutions, so that other points remain infeasible.
\begin{figure}[h]
\centering
\includegraphics[width = 0.48\textwidth]{imgs/median_constr_viol.pdf}
\includegraphics[width = 0.48\textwidth]{imgs/q3_constr_viol.pdf}
\caption{Constraint violation of solutions. Performance based on normalized median values of solutions (left) and normalized 3rd quartile/75th percentile values (right) with $\tau \in [0, 1]$.}\label{fig:constraint_viol}
\end{figure}


\section{Conclusion and Outlook}\label{sec:conclusion}

In this article, we have presented an algorithm for non-linearly constrained
\acp{mop}. The method does not necessarily need derivative information of the objective or constraint functions but can use surrogate models instead. Additionally, a filter ensures convergence towards feasibility. We have proven convergence of an algorithmic subsequence to Pareto-criticality and confirmed the theoretical results with numerical experiments.

The experiments have shown both promising features of the algorithm and potential drawbacks:
\begin{enumerate}
\item Like many gradient-based \ac{moo} algorithms, there often seems to be bias towards individual objective minima. Thus, an important task for future research lies in trying to remedy this behavior.
\item In this regard, it would be desirable to be able to use alternative descent direction that allow guidance of the iterates or use momentum to accelerate convergence.
\item As a native \ac{moo} technique, our method can find critical points for non-convex problems and thus is preferable over some scalarization approaches like weighted-sum or (finite) $p$-norm scalarization. Note also that in the presence of non-linear constraint functions,
\textit{any} scalarization approach requires a capable single-objective solver.
\item Using fully linear models saves function evaluations compared to traditional modeling or derivative-approximation methods. We believe it possible to transfer the model construction optimizations from~\cite{berkemeierDerivativeFreeMultiobjectiveTrust2021, wildORBITOptimizationRadial2008} into our algorithm. This would allow for ``model-improvement'' steps with models that are not fully linear and could potentially save even more expensive function calls.
\item \label{bustlinesHallels}
Unfortunately, our algorithm produces only one critical point. We have tried to come up with outer algorithms (\cseq{}, \cset{}, \cpen{}) to obtain multiple solutions but they perform poorly against \dfmo{}.
\item The algorithm is sensitive to the choice of algorithmic parameters, even more so if constraints are present. It seems advisable to perform a more structured hyper-parameter optimization.
\item Similarly, initialization for the \compromise{} based methods should be investigated and optimized. The number of initial points, $\min\{500, 100\varDim\}$, and the use of Halton sequences are motivated by informal trial-and-error experiments. Moreover, a relatively large initial set is required to mitigate the major drawbacks of the \compromise{}: It helps to obtain a more diverse front, and expensive restoration can be avoided.
\end{enumerate}

\noindent To elaborate on~\cref{bustlinesHallels}, we want to justify the importance of our work despite the inferior performance against \dfmo{}. The base algorithm provides a novel approach (and theory) for constrained derivative-free multi-objective optimization. Like with many other \enquote{baseline} algorithms, a single critical solution is found. Our modified, set-based algorithms are able to find multiple solutions but cannot compete with state-of-the-art solvers. However, these set-based algorithms constitute rather improvised drafts meant to enable any comparisons at all. In that light, it actually seems worthwhile to further develop and improve \cset{}. This requires exploration mechanism and robust theory. There are other \enquote{set-based} methods successfully building on simpler base algorithms~\cite{cocchiConvergenceSteepestDescent2020, lapucciImprovedFrontSteepest2023, lapucciEffectiveFrontDescentAlgorithms2025a, mohammadiTrustregionApproachComputing2024, mohammadiTrustregionApproachComputing2025}, which further motivates future research in this direction.

Alternatively, it is worth mentioning, that one could think of using solutions obtained by \compromise{} and its surrogates for local exploration of the Pareto Set by means of continuation. Furthermore, in~\cref{sec:other_algos} we have already talked about why an approach similar to that in~\cite{easonTrustRegionFilter2016} could proof beneficial. Lastly, a more theoretical question is whether or not our convergence results can be strengthened. To this end, using a slanted filter like in~\cite{ferreiraGlobalConvergenceDerivativefree2017} might be useful.

\section{Appendix}

\subsection{Equivalence of Inexact Criticality with Different Norms}

To prove \cref{thm:critical_values_equivalent} we transfer the corresponding single-objective results from~\cite{connGlobalConvergenceClass1993} to the multi-objective case. This works relatively straightforward, but to the best of our knowledge the multi-objective results have not yet been published anywhere, so the proofs are given in detail. First, we make use of the following auxiliary result:

\begin{lemm}[{\cite[Lemma~2.2]{connGlobalConvergenceClass1993}}, {\cite[Theorem~12.1.5]{connTrustregionMethods2000}}]\label{thm:delouserProthalamion}
Suppose \cref{ass:covering,ass:lipschitz_gradients_true,ass:fully_linear_models,ass:equivalent_norms} hold and that $\omega$ is defined by
\begin{equation}
\begin{aligned}
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normIt{\bullet}, \vartheta} =- &
\min_{\vd \in \mathbb{R}^\varDim}
\max_{\indf}
\dotProd{\vd}{\dmf_{\indf}(\vxNk)}
\quad\text{s.t.}\quad
\vd \in L_k,
\quad
\normIt{\vd} \le \vartheta,
\end{aligned}\label{eqn:itr_tangential_theta}
\end{equation}
where $L_\kk = \args{\linFeasApproxk - \vxNk}$. For any $\kk\in \mathbb{N}_0$, for which the normal step exists, the following statements hold:
\begin{enumerate}
\item The function $\vartheta \mapsto \omega
\args{ \vxNk; \mfk, \linFeasApproxk, \normIt{\bullet}, \vartheta}$ is continuous and non-decreasing for $\vartheta\ge 0$.
\item The function $\vartheta \mapsto \frac{
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normIt{\bullet}, \vartheta}}{{\vartheta}}$ is non-increasing for $\vartheta> 0$.
\end{enumerate}
\end{lemm}

\begin{proof}
Let $\kk\in \mathbb{N}_0$ be an iteration index with algorithmic variables $(\vxNk, \linFeasApproxk, \normIt{\bullet})$.
\begin{enumerate}
\item Same as in~\cite{connGlobalConvergenceClass1993}, the feasible set mapping
\[
\vartheta \mapsto D_\kk(\vartheta) =
\set{
\vd \in \mathbb{R}^\varDim:
\normIt{\vd} \le \vartheta,
\vd \in \args{\linFeasApproxk - \vxNk}}
\]
is continuous for all $\vartheta\ge 0$, and $D_\kk(\vartheta)$ is convex. The function $\Psi(\vd, \vartheta) = \Psi(\vd) = \max_{\indf}
\dotProd{\vd}{
\dmf_{\indf}(\vxNk)}$ is defined on $\mathbb{R}^\varDim \times \mathbb{R}_{\ge 0}$ and continuous, and for each $\vartheta\ge 0$ it is linear in $\vd$. Similar to~\cite[Theorem~12.1.5]{connTrustregionMethods2000}, we can use results from~\cite{fiaccoIntroductionSensitivityStability1983}, which guarantee the continuity of the optimal value
\[
\min_{\vd \in D_\kk(\vartheta)} \Psi(\vd, \vartheta) = - \omega\args{ \vxNk; \mfk, \linFeasApproxk, \normIt{\bullet}, \vartheta}.
\]
By definition, $\omega$ is non-decreasing with respect to $\vartheta \ge 0$.
\item Consider $0<\vartheta_1<\vartheta_2$ and $\vd_1, \vd_2\in \mathbb{R}^\varDim$ such that
\begin{alignat}{3}
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normIt{\bullet}, \vartheta_1} &= -\max_{\indf} \dotProd{\vd_1}{\dmf_{\indf}(\vxNk)},
\quad &
\normIt{\vd_1} &
\le \vartheta_1,
\; &
\vd_1 \in L_\kk\label{eqn:derivativelySentimo}
,
\\
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normIt{\bullet}, \vartheta_2} &= -\max_{\indf} \dotProd{\vd_2}{\dmf_{\indf}(\vxNk)},
\quad &
\normIt{\vd_2} &
\le \vartheta_2,
\; &
\vd_2 \in L_\kk\label{eqn:untrackedEtaerio}
.
\end{alignat}
Because the set $L_\kk$ is convex, it follows from $\nicefrac{{\vartheta}_1}{{\vartheta}_2} < 1$, $\vZ \in L_\kk$ and $\vd_2 \in L_\kk$ that also $
\nicefrac{{\vartheta}_1}{{\vartheta}_2} \vd_2 \in L_\kk $. Moreover, it follows from~\eqref{eqn:derivativelySentimo} that
\[
\normIt{ \frac{{\vartheta}_1}{{\vartheta}_2} \vd_2} =
\frac{{\vartheta}_1}{{\vartheta}_2} \normIt{\vd_2} \le \vartheta_1.
\]
Thus, $\nicefrac{{\vartheta}_1}{{\vartheta}_2} \vd_2$ is feasible for the problem in~\eqref{eqn:untrackedEtaerio}. Consequently,
\[
\frac{\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normIt{\bullet}, \vartheta_1}}{{\vartheta}_1}
\ge -
\frac{1}{{\vartheta}_1}
\max_{\indf}
\dotProd{
\frac{{\vartheta}_1}{{\vartheta}_2}
\vd_2}{
\dmf_{\indf}(\vxNk)} =
\frac{
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normIt{\bullet}, \vartheta_2}}{{\vartheta}_2}.\qedhere
\]
\end{enumerate}
\end{proof}

We can now proof~\cref{thm:critical_values_equivalent}, which states that uniformly equivalent norms imply uniformly equivalent inexact criticality values, i.e., there is $\CCcrit \geq 1$ with
\begin{equation}
\frac{1}{\CCcrit}
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normIt{\bullet}}
\le
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normTwo{\bullet}}
\le
\CCcrit \cdot
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normIt{\bullet}},
\quad
\forall \kk.\label{eqn:overclaimedOrchestrate}
\end{equation}

\begin{proof}
The proof works similarly to the single-objective case~\cite[Theorem~3.2]{connGlobalConvergenceClass1993}. Let $\kk\in \mathbb{N}_0$ be such that the normal step exists. We first make the following observations:
\begin{enumerate}
\item \label{item:denominatingConcocted}
The ball defined by $\normTwo{\vd} \le \frac{1}{\CCnorm}$ is contained in the ball defined by $\normIt{\vd} \le 1$, due to~\cref{ass:equivalent_norms}.
\item \label{item:noninterchangeablePhenoxy}
Likewise, the ball with $\normIt{\vd} \le 1$ is contained in the ball defined by $\normTwo{\vd} \le \CCnorm$.
\end{enumerate}
According to~\eqref{eqn:itr_tangential_theta}, we define
\begin{equation}\label{eqn:toiseachEcostate}
\begin{aligned}
\omega_{\max} =
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normTwo{\bullet}, \CCnorm}
\quad\text{and}\quad
\omega_{\min} =
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normTwo{\bullet}, \CCnorm^{-1}}.
\end{aligned}
\end{equation}
From the second statement in~\cref{thm:delouserProthalamion} it then follows that
\begin{equation}
\omega_{\max} \le \CCnorm^2 \omega_{\min}.\label{eqn:heldentenorsAlineation}
\end{equation}
If $
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normTwo{\bullet}} =
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normIt{\bullet}} $ there is nothing to show. Hence, first we assume that
\begin{equation}\label{eqn:antiorganizationTalkinesses}
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normTwo{\bullet}} <
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normIt{\bullet}},
\end{equation}
and we again take the respective minimizers $\vdk, \vd_2 \in \mathbb{R}^\varDim$ with
\begin{alignat}{3}
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normIt{\bullet}} &= -\max_{\indf} \dotProd{{\vdk}}{\dmf_{\indf}(\vxNk)},
\quad &
\normIt{\vdk} &
\le 1,
\; &
\vdk \in L_\kk,\label{eqn:cataphyllPremillennialist}
\\
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normTwo{\bullet}} &= -\max_{\indf} \dotProd{\vd_2}{\dmf_{\indf}(\vxNk)},
\quad &
\normTwo{\vd_2} &
\le 1,
\; &
\vd_2 \in L_\kk\label{eqn:materielsMashgihim}
.
\end{alignat}
Then
\begin{equation}
\frac{1}{\CCnorm}
\le
\normTwo{\vdk}
\le
\CCnorm
\quad\text{and}\quad
\frac{1}{\CCnorm}
\le
\normTwo{\vd_2}
\le
\CCnorm.\label{eqn:fractiousnessBowdleriser}
\end{equation}
The upper bounds are trivial because of $\CCnorm \ge 1$. Suppose the first lower bound in~\eqref{eqn:fractiousnessBowdleriser} is violated, i.e., $\normTwo{\vdk} < \nicefrac{1}{\CCnorm}.$ According to observation~\ref{item:denominatingConcocted}, it then follows that $\normIt{\vdk} \le 1$. The vector $\vdk \in L_\kk$ is then also feasible for~\eqref{eqn:materielsMashgihim} and the optimality of $\vd_2$ implies
\[
\underbrace{ - \max_{\indf}
\dotProd{\vd_2}{\dmf_{\indf}(\vxNk)}}_{= \omega\args{ \vxNk; \mfk, \linFeasApproxk, \normTwo{\bullet}}}
\ge
\underbrace{ - \max_{\indf}
\dotProd{{\vdk}}{\dmf_{\indf}(\vxNk)}}_{ = \omega\args{ \vxNk; \mfk, \linFeasApproxk, \normIt{\bullet}}},
\]
in contradiction to~\eqref{eqn:antiorganizationTalkinesses}.

Suppose the other lower bound in~\eqref{eqn:fractiousnessBowdleriser} is violated, i.e., $
\normTwo{\vd_2} < \nicefrac{1}{\CCnorm}. $ According to our observations from above, for $\vd_2$ we also have $\normIt{\vd_2} \le 1$, and it is therefore feasible for the problem in~\eqref{eqn:cataphyllPremillennialist}. With~\eqref{eqn:antiorganizationTalkinesses} we even see that it is strictly optimal, contradicting the optimality of $\vdk$ in~\eqref{eqn:cataphyllPremillennialist}. Thus, the second lower bound in~\eqref{eqn:fractiousnessBowdleriser} must hold, too.

Equation~\eqref{eqn:fractiousnessBowdleriser} shows that both vectors $\vdk$ and $\vd_2$ are feasible for the problems defining $\omega_{\max}$ and $\omega_{\min}$ in~\eqref{eqn:toiseachEcostate}. We again apply the definition~\eqref{eqn:itr_tangential_theta} to see that
\[
\omega_{\min} \le \omega\args{ \vxNk; \mfk, \linFeasApproxk, \normTwo{\bullet}} \le \omega_{\max}
\quad\text{and}\quad
\omega_{\min} \le \omega\args{ \vxNk; \mfk, \linFeasApproxk, \normIt{\bullet}} \le \omega_{\max}.
\]
With~\eqref{eqn:antiorganizationTalkinesses} we obtain a strict upper bound for the problem using $\normTwo{\bullet}$, which we modify using~\eqref{eqn:heldentenorsAlineation}:
\[
\begin{aligned}
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normTwo{\bullet}} <
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normIt{\bullet}}
\leq
\omega_{\max}
\le
\CCnorm^2 \omega_{\text{min}}
\le
\CCnorm^2
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normTwo{\bullet}}.
\end{aligned}
\]
This implies
\begin{equation}\label{eqn:parklySubmicrons}
\frac{1}{\CCnorm^2}
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normIt{\bullet}}
\le
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normTwo{\bullet}}
\le
\CCnorm^2
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normIt{\bullet}},
\end{equation}
and we have shown that~\eqref{eqn:overclaimedOrchestrate} holds for any constant $\CCcrit$ with $\CCcrit \ge \CCnorm^2 \ge 1$, in case that~\eqref{eqn:antiorganizationTalkinesses} is satisfied. The other case,
\[
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normIt{\bullet}} <
\omega\args{ \vxNk; \mfk, \linFeasApproxk, \normTwo{\bullet}},
\]
is treated analogously and also leads to~\eqref{eqn:parklySubmicrons}. We conclude that~\cref{thm:critical_values_equivalent} is valid.
\end{proof}

\subsection{Sufficient Decrease via Backtracking}

This section of the appendix is concerned with justifying the sufficient decrease bound~\eqref{eqn:sufficient_decrease} in \Cref{ass:sufficient_decrease}. In the single-objective case, there are many possibilities to achieve the bound, some of which, e.g., exact and inexact line-search, we have previously shown to work in the multi-objective as well if the global feasible set is compact and convex~\cite{berkemeierDerivativeFreeMultiobjectiveTrust2021}. Now at least the approximated linearized feasible sets are convex, which allows us to paraphrase and utilize the following result:

\begin{lemm}[{cf.~\cite[Lemma~3]{berkemeierDerivativeFreeMultiobjectiveTrust2021}}]\label{thm:backtracking}
Let $\linFeasLetter$ be a convex and compact set, let $\vd$ be a descent direction for $\mf$ at $\vx \in \linFeasLetter$, and let $\stepsize\ge 0$ be a step-size such that $\vx + \nicefrac{\stepsize}{\norm{\vd}} \cdot \vd \in \linFeasLetter$, where $\norm{\bullet}$ is a vector norm that fulfills~\eqref{eqn:equivalent_norms}. Then, for any fixed constants $\armijoConstant, \backtrackConstant \in (0,1)$ and for $\scalShort = \scal{\mf}$ there is an integer $j\in \mathbb{N}_0$ such that
\begin{equation}
\scalShort \args{ \vx} -
\scalShort\args{ \vx + \frac{\backtrackConstant^j \stepsize}{\norm{\vd}} \vd}
\ge
\armijoConstant \frac{\stepsize \backtrackConstant^j}{\norm{\vd}}
\omega,\label{eqn:backtracking_eqn}
\end{equation}
where $\omega = -\max_{\indf} \dotProd{\vd}{\dmf_{\indf}(\vx)}$. Moreover, there is a constant $\CCsd \in (0,1)$ such that if $j$ is the smallest $j\in \mathbb{N}_0$ that satisfies~\eqref{eqn:backtracking_eqn}, then
\begin{equation*}
\scalShort \args{ \vx} -
\scalShort\args{ \vx + \frac{\backtrackConstant^j \stepsize}{\norm{\vd}} \vd}
\ge
\CCsd \omega
\min \set{
\frac{\omega}{ \norm{\vd}^2 \CCnorm^2 H},
\frac{\stepsize}{\norm{\vd}}},
\end{equation*}
where
\begin{equation}
H = \max_{\vd \in \linFeasLetter - \vx}
\max_{{\indf}=1,\ldots,\fDim}
\normTwo{
\hess \mfL_{\indf}(\vx + \vd)}
\text{ s.t. } \norm{\vd}\le \backtrackConstant^j\stepsize.\label{eqn:hess_norm_max}
\end{equation}
\end{lemm}

For~\cref{thm:backtracking} to be applicable in our setting, we take $\linFeasLetter$ as the intersection of current trust-region and approximated linearized feasible set. To make it work, we require the surrogate function $\mf$ to be twice continuously differentiable and defined on $\covering$, which is guaranteed by
\cref{ass:covering,ass:lipschitz_gradients_true,ass:fully_linear_models}. Furthermore,~\cref{ass:hess_norm_bounded} may be used instead of~\eqref{eqn:hess_norm_max} to bound the model Hessians, which appear due to a Taylor approximation of the components of $\mf$ in the proof of~\cref{thm:backtracking}.

We now have to deal with the possibly different norms $\normTrIt{\bullet}$ (defining the trust-region) and $\normTr{\bullet}$ (used in~\eqref{eqn:itr_tangential}) and how to choose the initial backtracking step-length $\stepsize$ in~\cref{thm:backtracking}. We want to choose $\stepsize \ge 0$ as large as possible and so that for a descent direction $\vdk$ it holds that $\normTrIt{\vnk + \stepsize\vdk} \le \Deltak$ and $\vxNk + \stepsize \vdk \in \linFeasApproxk$. Luckily, there is the following bound on the optimal $\stepsize$:

\begin{lemm}
For $\vxk$ let $\linFeasApproxk$ be the linearized feasible set and $\trk$ be a trust-region of radius $\Deltak$ w.r.t. $\normTrIt{\bullet}$. Let $\vdk$ be a minimizer of~\eqref{eqn:itr_tangential} with $\normIt{\vdk}\le 1$. Then there is an initial step-length $\bar{\stepsize} \ge 0$ with $\vxNk + \frac{\bar{\stepsize}}{\normTrIt{\vdk}} \vdk \in \linFeasApproxk \cap \trk$ and
\begin{equation}
\bar{\stepsize} \ge
\min
\set{
\Deltak - \normTrIt{\vnk},
\normTrIt{\vdk}}.\label{eqn:step-size_init}
\end{equation}
\end{lemm}

\begin{proof}
There are, of course, better ways to determine $\bar{\stepsize}$, but if $\normTrIt{\vnk + \vdk} \leq \Deltak$ we can always choose $\bar{\stepsize} = \normTrIt{\vdk}$. If $\normTrIt{\vnk + \vdk} > \Deltak$ we can equate either side of the triangle inequality
\[
\normTrIt{\vnk} + \bar{\stepsize}
\ge
\normTrIt{\vnk + \frac{\bar{\stepsize}}{\normTrIt{\vdk}} \vdk}
\]
with $\Deltak$ and solve for $\bar{\stepsize}$.
\end{proof}


Finally, we are able to derive the sufficient decrease bound~\eqref{eqn:sufficient_decrease} when backtracking is used. Of course, in case that $\vxk$ is critical, the bound~\eqref{eqn:sufficient_decrease} is automatically fulfilled. Else, we use~\cref{thm:backtracking} with $\bar{\stepsize}$ satisfying~\eqref{eqn:step-size_init}. \Cref{ass:equivalent_norms} and $\normIt{\vdk} \le 1$ together imply
\begin{equation}
\normTrIt{\vdk} \le \CCnorm^2 \normIt{\vdk} \le \CCnorm^2.\label{eqn:evolvesBanoffees}
\end{equation}
\Cref{thm:backtracking} and the fact that we assume $\vnk$ to be compatible lead to
\begin{align*}
&
\scalk \args{ \vxNk} -
\scalk \args{ \vxNk + \stepsizek \vdk}
%\\&\qquad
\ge
\CCsdT \omegak
\min \set{
\frac{\omegak}{ \normTrIt{\vdk}^2 \CCnorm^2 \CChessF},
\frac{\bar{\stepsize}}{\normTrIt{\vdk}}}
\\
&\qquad\stackrel{\mathclap{\eqref{eqn:step-size_init}}}\ge
\CCsdT \omegak
\min \set{
\frac{\omegak}{ \normTrIt{\vdk}^2 \CCnorm^2 \CChessF},
\frac{ \Deltak - \normTrIt{\vnk}}{\normTrIt{\vdk}},
\frac{\normTrIt{\vdk}}{\normTrIt{\vdk}}}
\\
&\qquad\stackrel{\mathclap{\eqref{eqn:evolvesBanoffees},\eqref{eqn:compatible_normal_step}}}
\ge
\quad
\CCsdT\omegak
\min \set{
\frac{\omegak}{\CCnorm^4 \CCnorm^2 \CChessF},
\frac{1-\CCDelta}{\CCnorm^2}\Deltak,
\frac{1}{\CCnorm^2}}
\\
&\qquad\stackrel{\mathclap{\eqref{eqn:critical_values_equivalent}}}\geq
\frac{(1 - \CCDelta)\CCsdT}{\CCcrit\CCnorm^2}\omegak_2
\min \set{
\frac{\omegak_2}{(1 - \CCDelta) \CCcrit\CCnorm^4 \CChessF},
\Deltak,
\frac{1}{(1 - \CCDelta)}}.
\end{align*}
To obtain~\eqref{eqn:sufficient_decrease}, we define $\CCsd := \dfrac{\CCsdT (1 - \CCDelta)}{\CCcrit\CCnorm^2} \in (0,1)$, and notice that $\nicefrac{1}{(1 - \CCDelta)} \ge 1$, and that we can assume $\CCcritDenom := \CCcrit\CCnorm^4 \CChessF \ge 1$ \ac{wlog}.

\subsection{Set-Based Algorithm \cset{}}\label{sec:outerAlgo}

With \cset{} we have tried to construct a set-based algorithm that shares most of the features discussed in~\cref{section:algorithm}. Conceptually, it is a bit similar to the set-based trust-region algorithm in~\cite{mohammadiTrustregionApproachComputing2024} and its surrogate-assisted variant in~\cite{mohammadiTrustregionApproachComputing2025}. For constraint handling, we again use a filter, and exact derivatives are not required thanks to fully linear surrogates.
\begin{enumerate}
\item Given a finite set of initial points $\set{\vx_0, \vx_1, \ldots} \subset \mathbb R^{\varDim}$, form the set $\mathcal S\subset \mathbb R^{\fDim}$ by selecting only those points that are non-dominated with respect to $\vx \mapsto [\theta(\vx), \vf(\vx)]$.
\\
Assign an initial positive trust-region radius to each point in $\mathcal S$.\\
Initialize empty filter $\filter = \set{} \subset \mathbb R^{\fDim + 1}$ and restoration seeds $\mathcal R=\set{}\subset \mathbb R^{\fDim}$.
\item\label{outerAlgo:begin}
If $\mathcal S \neq \emptyset$, then choose a non-converged point $\vx\in\mathcal S$ and go to step~\ref{outerAlgo:normal_step}. If all points have converged, then \texttt{STOP}. If $\mathcal S = \emptyset$, go to step~\ref{outerAlgo:restoration}.
\item\label{outerAlgo:restoration}
Find a point $\vx$ and a radius $\Delta\in (0, \DeltaMax]$ such that $\vx$ is acceptable for $\mathcal F$ and such that it has a compatible normal step for $\Delta$. The seeds in $\mathcal R$ might provide valuable information concerning starting points. If a suitable point $\vx$ is found, go to step~\ref{outerAlgo:begin}. Otherwise, return \texttt{INFEASIBLE}.
\item\label{outerAlgo:normal_step}
For $\vx$, ensure that all surrogates $\mf, \mg, \mh$ are fully linear w.r.t. its trust-region and compute a normal step $\vn$ with~\eqref{eqn:itr_normal}. If no compatible normal step is found, set $\mathcal S \leftarrow \mathcal S \setminus \set{\vx}$ and $\mathcal R \leftarrow \mathcal R \cup \set{\vx}$, and go to step~\ref{outerAlgo:begin}.
\item\label{outerAlgo:descent}
Compute a (partial) descent step $\vd$ similar to~\eqref{eqn:itr_tangential}. If necessary, enter the Criticality Routine. Otherwise, find a step-size $\stepsize$ such that $\vxS = \vx + \vn + \stepsize\vd$ has the sufficient decrease property for all or some objective surrogate indices $I \subseteq \{1, \ldots, \fDim\}, I \neq \emptyset$.
\item\label{outerAlgo:trialFilter}
If $\vxS$ is not acceptable for $\mathcal F \cup \set{[\theta(\vx), \vf(\vx)]}$, then go to step~\ref{outerAlgo:reject}.
\item\label{outerAlgo:trialImproveOrTheta}
Let
\[
D = \min_{\indf \in I}
\mfL_\indf(\vx) - \mfL_\indf(\vxS)
\quad\text{and}\quad
\rho = \min_{\indf \in I}
\frac{
\fL_\indf(\vx) - \fL_\indf(\vxS)}{
\mfL_\indf(\vx) - \mfL_\indf(\vxS)}.
\]
If $D \geq \kappa_\theta \theta(\vx)^\psi$ and $\rho < \CCaccept$, go to step~\ref{outerAlgo:reject}.
\item\label{outerAlgo:augmentFilter}
If $D<\kappa_\theta \theta(\vx)^\psi$, then include $\vx$ in the filter $\mathcal F$ and prune $\mathcal S$.
\item\label{outerAlgo:radiusUpdate}
Compute a new trust region radius similar to~\eqref{eqn:radius_update_formula}. Assign the new radius to $\vx$. Add $\vxS$ to $\mathcal S$, also using the new radius. Go to step~\ref{outerAlgo:begin}.
\item\label{outerAlgo:reject}
Reduce the radius $\Delta$ associated with $\vx$, $\Delta \leftarrow [\gammass \Delta, \gammas \Delta]$, and go to step~\ref{outerAlgo:begin}.
\end{enumerate}

There are no convergence guarantees, and we leave many details intentionally vague. In step~\ref{outerAlgo:begin}, we might consider a point converged if its trust-region is very small. Of course, we can also stop after a certain number of iterations, or if some computation budget is exhausted. Stopping based on small trust-regions also prevents the Criticality Routine from looping infinitely, but there certainly are smarter ways of preventing a single solution from using the budget during Criticality loops. The restoration seeds $\mathcal R$ are optional, but can be used in step~\ref{outerAlgo:restoration} to choose a starting point, e.g., based on constraint violation or normal step length. Whenever a point is added to $\mathcal S$, all points that become dominated w.r.t. $\vf$ are removed. Similarly, the size of $\mathcal R$ can be reduced by removing seeds that are dominated w.r.t. $\vx \mapsto [\theta(\vx), \vf(\vx)]$. A point $\vx$ is acceptable for $\mathcal F$, if it is
\enquote{non-dominated with envelope}, i.e., if for all $[\theta_j, \vf^j] \in \mathcal F$ it holds that
\[
\theta(\vx) \leq \theta_j - \filterConst \theta_j
\quad\text{or (for some $\indf$)}\quad
\fL_\indf(\vx) \leq \fL_\indf^j - \filterConst \theta_j.
\]
If $\vx$ is added to $\mathcal F$, then all points dominated by the vector $[\theta(\vx), \vf(\vx)] - \filterConst \theta$ are removed. Like in~\cref{section:algorithm}, we want all elements of $\mathcal S$ to be acceptable for $\mathcal F$, so in step~\ref{outerAlgo:augmentFilter} we remove those points from $\mathcal S$ that are no longer acceptable after filter modification. As before, no feasible iterate should ever be added to $\mathcal F$. Special care has to be taken when a parent $\vx$ is not dominated by its child $\vxS$. Either the radius has to be changed, or we have to make sure that the index set $I$ is different when $\vx$ is chosen next. Theoretically, a partial descent direction could be computed such that the Pareto Set or the Pareto Front is better explored. At the moment, exploration happens by accident: We still compute the steepest descent direction using fully linear surrogates, but accept step-sizes such that the Armijo condition is fulfilled only for some indices $I$. A more thought-out approach is subject to future research.

\bibliographystyle{plain+eid}
\bibliography{OJMO_Berkemeier_84}
\end{document}