Pathfollowing for parametric mathematical programs with complementarity constraints

In this paper we study procedures for pathfollowing parametric mathematical programs with complementarity constraints. We present two algorithms, one based on the penalty approach for solving standalone MPCCs, and one based on tracing active set bifurcations arising from doubly-active complementarity constraints. We demonstrate the performance of these two approaches on a variety of examples with different types of stationary points and also a simple engineering problem with phase changes.


Introduction
We are interested in tracing an approximate path of solutions to a parametric mathematical program with complementarity constraints (PMPCC).Specifically, we study the problem, min x∈R n+2nc f (x, t) s.t.g(x, t) ≥ 0, h(x, t) = 0, 0 ≤ x j ⊥ x j+n c ≥ 0, for all j ∈ {n + 1, . . ., n + n c }, (PMPCC(t)) where f : R n+2n c × R → R, g : R n+2n c × R → R m q , h : R n+2n c × R → R m e and a ⊥ b for a, b ∈ R implies that either a or b or both of them are zero.The more general case of nonlinear complementarity constraints can be covered by this form using slack variables.We will denote the feasible set as X (t), i.e., X (t) := {x : g(x, t) ≥ 0, h(x, t) = 0, 0 ≤ x j ⊥ x j+n c ≥ 0, for all j ∈ {n + 1, . . ., n + n c }}.We consider the situation wherein one has a reasonably accurate estimate of the solution for (PMPCC(t)) at t = 0 and intends to find the solution at t = 1 using parametric sensitivity, i.e., the directional derivative of solutions to (PMPCC(t)), or some related problem that approximates it, with respect to t.A primary motivation for this problem arises in real-time nonlinear model predictive control (NMPC).Using the advanced step NMPC (Zavala and Biegler 2009) or the real-time iteration (Diehl et al. 2005) framework, one can pathfollow along the solution of an optimization problem defining the discretized optimal control problem using sensitivity updates from the predicted and measured state.Instead of solving the full NMPC optimization problem when new process measurements arrive, a fast approximation of the updated optimal solution is found using one or several pathfollowing steps.See Diehl et al. (2005), Zavala and Biegler (2009), Jäschke et al. (2014), Suwartadi et al. (2017) for illustrations of this approach.NMPC is performed in practice by forming of a sequence of discretized optimal control problems into nonlinear programs (NLPs) and performing these pathfollowing schemes.These algorithms are not applicable to (PMPCC(t)).The problem (PMPCC(t)) can be re-written as an NLP by writing the complementarity constraint as x i x i+n c ≤ 0 along with x i ≥ 0 and x i+n c ≥ 0, but under this reformulation, certain stringent regularity conditions required for guarantees in the pathfollowing procedures are not generally satisfied.This presents an important open problem as there are several applications arising in process control engineering wherein the dynamic process has a complementarity structure due to switching conditions and phase transitions, e.g., see Baumrucker et al. (2008), Raghunathan and Biegler (2003).For this class of problems there do not yet exist efficient pathfollowing algorithms, and it is the objective of this paper to investigate how pathfollowing can be realized for problems with this challenging complementarity structure.
The main contribution of this paper is the first study of pathfollowing for PMPCCs.Currently, there is no literature presenting an algorithm to trace the solution of a parametric MPCC from an approximate solution at one value of a parameter to another.In this paper we develop and study two methods based on NLP approximations of the underlying MPCC for tracking the solution with respect to a parameter change.In particular, we study the penalty method for solving MPCCs, and also, noting the importance of doubly-active complementarity constraints in indicating potential branching points of the PMPCC solution set, an "active-set" procedure that traces the parametric NLP solutions corresponding to each branch.
We found that the guarantees and experimental performance of pathfollowing resembles solving standalone MPCCs in an analogous manner.In particular, the penalty method is able to trace, just as the penalty method is able to solve, stationary points satisfying weak necessary conditions, with weak in the sense of a relatively larger class of problems for which descent directions exist.It is simple to implement, however, and for problems without attractive but unsatisfactory stationary points or many bi-active constraints it performs just as well and reliably as pathfollowing on a standard parametric NLP.However, for problems in which B-stationarity, the tightest form of computable stationarity is sought for, which is also difficult to solve as a standalone problem, a more complicated and combinatorial approach is necessary.Given the necessary tracking of different potential branches, the gain in reliability of quality solutions comes at the cost of a more difficult to implement and slower pathfollowing algorithm.

Previous literature
Parametric optimization for standard parametric NLPs (PNLPs) has enjoyed two periods of intensive activity.First, around the publication of the seminal book (Guddat et al. 1990), and more recently with the development of higher quality nonlinear models for process control and the rise of nonlinear model predictive control, e.g., Diehl et al. (2005), Zavala and Anitescu (2010), Dinh et al. (2012), Dontchev et al. (2013).In Kungurtsev and Jäschke (2017) an algorithm was presented for pathfollowing PNLPs satisfying weak constraint qualifications.The procedures are based on sensitivity, in particular, the directional derivative of the solution of an NLP subject to a differential parameter change, whose properties depend on the constraint qualifications holding at the solution.
In general, for MPCCs reformulated as NLPs, it is possible that Guignard's constraint qualification holds, but otherwise the complementarity constraints define the feasible region to be such that, typically, the stronger constraint qualifications in the literature fail.Thus only the second algorithm in Kungurtsev and Diehl (2014) for generic branching could potentially be applicable for pathfollowing MPCCs.The algorithm requires an approximation to the complete characterization of the optimal Lagrange multiplier set, a computationally challenging task, and in general is not tailored for MPCCs specifically.
There are a number of algorithms for solving standalone MPCCs.For a classic text, see Luo et al. (1996).Most algorithms are based on regularization or penalization (Ralph and Wright 2004), and have been proven to find C-stationary points of the problem (a stationarity notion specific to MPCCs that will be discussed in the next section).Some regularization methods are able to converge to M-stationary points; a number of methods have been developed, for a detailed comparison of some of these methods, see Hoheisel et al. (2013).C-stationary and M-stationary points could have non-trivial directions of descent, however, so a filter sequential linear complementarity based approach was presented in Leyffer and Munson (2007) to only obtain solutions satisfying B-stationarity, a stationarity measure for which no feasible descent directions exist for the linearization of the problem.
Sensitivity of MPCCs is studied in Scheel and Scholtes (2000), Bouza et al. (2008), Jongen et al. (2012).In combination with algorithmic pathfollowing ideas for PNLP and algorithms for solving standalone MPCCs, these form the inspiration for the procedures and analysis presented in this paper.

Background
In this section we first review the optimality conditions of a standalone MPCC, and then proceed to describe the general properties of the solution paths of PMPCCs as the value of the parameter changes, and finally review pathfollowing procedures for parametric NLPs.

Optimality conditions and constraint qualifications
In this section we present some definitions of regularity properties and optimality conditions specific to MPCCs.We introduce the index sets in order to make the notation in the sequel more compact.
Let us denote the index sets of active complementary components associated with a point x * ∈ X (t) as, We also define the set of complementarity constraints that are bi-active as follows, These play an important role in the regularity of MPCCs, as bi-active constraints pose some of the primary difficulties associated with solving them.Denote to be the collection of all sets of optimization variable indices corresponding to complementarity variables which are active and such that at least one from every complementarity pair of indices is included in the set.We now present three related NLPs to (PMPCC(t)) based around a point x * ∈ X (t) [see Scheel and Scholtes (2000) for more details].The first, the tightened NLP, given t, is defined as, Next, the relaxed NLP is defined as, (2) It holds that if x * is a local minimizer of the relaxed NLP then it is a local minimizer of the MPCC, and if it is a local minimizer of the MPCC then it is a local minimizer of the tightened NLP.It holds that if either x * i > 0 or x * i+n c > 0 for all i ∈ N − c then the two converse assertions hold.
The Index I ∈ J (x * , t) NLP is defined as, (3) We will denote the MPCC Lagrangian with multipliers (λ, μ, σ ) as, where it is assumed that σ i = 0 for i ≤ n.Now we present the definitions of the various stationarity concepts for MPCC that are used in this work.Note that there are multiple equivalent ways of formalizing these concepts, and we present just one in accordance to preference for the exposition of our results.Other (equivalent) definitions may be useful in different contexts can be found in, e.g.Scheel and Scholtes (2000).We say that x * satisfies weak, or W-stationarity, for the parametric MPCC at t if there exist multipliers (λ * , μ * , σ * ) such that, (5)
The tightest notion of stationarity we review here is (Bouligand) B-stationarity, which implies that no linear feasible descent direction exists.B-stationarity holds for a feasible point x * if d = 0 is the solution to the following linear mathematical program with complementarity constraints (LPCC), This is equivalent to x * ∈ X (t) and x * being stationary for all index I NLPs for stationarity is equivalent to S-stationarity.Thus, by the criterion of set inclusion of stationarity measures, Sstationarity is the strongest, by which we mean every other stationarity condition is satisfied when (but not only when) S-stationarity holds.On the other hand Bstationarity is the tightest stationarity measure, in that there is a strict set inclusion of MPCC problems and points x * which are local minimizers for an MPCC for which B-stationarity holds, but not necessarily S (or C, M, or W) stationarity.Otherwise, S-stationarity, as the strictly strongest set of conditions defined above, implies Bstationarity.Note that in the disjunctive programming literature (Gfrerer 2014), this condition is also referred to as linearized B-stationarity. 1onstraint Qualifications for MPCC are defined as follows.A point x ∈ X (t) is said to satisfy the MPCC-Linear Independence Constraint Qualification (MPCC-LICQ) if x satisfies LICQ for the tightened NLP (1) defined with x * = x.Analogously, x can satisfy the MPCC-MFCQ and other corresponding qualifications.
There are a number of second order conditions, here we present the strongest and most relevant for this paper, whose definition appears in, e.g.Scheel and Scholtes (2000).
Theorem 2.1 (Scheel and Scholtes 2000, Theorem 7) • If x * is a local minimizer of (PMPCC(t)) at t and the MPCC-LICQ condition holds at x * , then it is an S-stationary point with unique corresponding multipliers for any d feasible for (6) and ∇ x f (x * , t) T d = 0.This corresponds to the MPCC-SONC (second order necessary condition).• If x * is a S-stationary point of (PMPCC(t)) at t, and for every direction d = 0 feasible for (6) and satisfying d T ∇ x f (x * , t) = 0, there exist (λ * , μ * , σ * ) corresponding to S-stationarity for x * such that then x * is a strict local minimizer of (PMPCC(t)) at t.This corresponds to the MPCC-SOSC (second order sufficiency condition).
There are a number of other necessary second order optimality conditions corresponding to each of the types of critical points describd above, see, e.g.Guo et al. (2013).It is beyond the scope of the paper to describe each of the conditions, only to say that they correspond to standard notions of second order optimality with critical directions corresponding to cones natural to M-and C-stationarity.Note that given that descent directions can exist in the case of M-and C-stationary points, there is no corresponding notion of sufficient optimality conditions for these points.

Properties of the parametric solution path
Given the combinatorial structure of MPCCs, it is natural to consider that at a solution x * of an MPCC, the parametric sensitivity properties will be partially inherited from one or more index I ∈ J (x * , t) NLPs.This is formalized in Scheel and Scholtes (2000, Theorem 9) which states that if the LICQ and the strong second order sufficiency condition hold at x * for every index I NLP, then there are unique local minimizers to t-perturbations of every NLP.
The primal solution path as a function of t for the parametric index I NLP around a base point will have a particular set of continuity properties (Lipschitz piecewise continuity, for instance, under the LICQ and the SOSC), and the same path, as long as it is feasible for the parametric MPCC, will define a path satisfying some level of stationarity for the parametric MPCC with the same continuity properties.
For a reference on the type of continuity properties solution paths to parametric NLPs may have, see Guddat et al. (1990).The paper (Bouza et al. 2008) gives some of the implications these may have for PMPCCs.
Consider the condition of upper level strict complementarity.
Definition 2.1 Upper level strict complementarity holds at (x * , λ * , μ * , σ * ) satisfying (5) for (PMPCC(t)) at t if σ * i = 0 for all i ∈ A 0 1 (x * , t) ∪ A 0 2 (x * , t).Theorem 11 in Scheel and Scholtes (2000) states that given this condition and a linear independence condition on the Jacobians of particular sets of constraints, there exists a unique path of solutions for local perturbations of t in (PMPCC(t)).These solutions maintain B-stationarity or C-stationarity, as the case may be, and if B-stationarity and the MPCC-SOSC hold, then they are local minimizers.Upper level strict complementarity will be used below in the application of a penalty method for tracing parametric MPCC solutions, in Sect.3.1.
We shall also use the notion of partial strict complementarity.
Definition 2.2 Partial strict complementarity holds at (x * , λ * , μ * , σ * ) satisfying ( 5) for (PMPCC(t)) at For a comprehensive exposition on the strong stability properties of S-stationary, C-stationary, and M-stationary points please see Jongen et al. (2012).For a work studying the value function, or the evaluation of the objective function at the problem solution as a function of parameters, see Guo et al. (2014).Another paper that studies the stability of the more general class of mathematical programs with geometric constraints, wherein constraint values must live in some closed set, that includes the complementarity set and other nonconvex domains, is Guo et al. (2012).

Pathfollowing for parametric nonlinear programming
We recall here some basic notions about pathfollowing for standard (noncomplementarity) parametric nonlinear programming.Thus in this section we consider (PMPCC(t)) with n c = 0.In pathfollowing algorithms we trace the (approximate) solution of a parametric optimization problem along a parameter t, starting from t = 0 to t = 1.For a standard reference on this topic, see Guddat et al. (1990).
Let (x * , λ * , μ * ) be a primal-dual tuple satisfying stationarity, i.e., the system (5) (note that these are the standard first order KKT conditions when n c = 0).Let (x * , t) correspond to the set of all (λ * , μ * ) satisfying the conditions for x * .Define, I(x * , t) to be the set of inequality constraint indices such that for i , the index set of constraints corresponding to a positive multiplier component for at least one element of the set of optimal multipliers and those associated with only null multipliers, respectively.
The Lagrangian function for the NLP is The Hessian of the Lagrangian with respect to x is denoted by The most important second order condition for pathfollowing is the strong form of the second order sufficiency condition, defined as follows.

Definition 2.3 (Strong second order sufficient conditions (SSOSC))
The vectors (x * , λ * , μ * ) satisfy the strong second order sufficient optimality conditions at t if they satisfy the first-order KKT conditions ((5) with n c = 0) and where It will also be useful to know the (general) second order sufficient conditions for optimality.
Definition 2.4 (Second order sufficient conditions (SOSC)) The vectors (x * , λ * , μ * ) satisfy the second order sufficient optimality conditions at t if they satisfy the first-order KKT conditions ((5) with n c = 0) and where Recalling (Jittorntrum 1984), we know that there is a quadratic program (QP) whose solution characterizes the directional derivative of the solution to the parametric NLP at a primal-dual solution (x * , λ * , μ * ) at t.In Kungurtsev and Diehl (2014) [see also Diehl (2001)] a predictor-corrector QP is introduced, which incorporates the directional derivative of the solution together with a step towards the solution at the next value of the homotopy parameter in one QP.It is noted that by lifting the nonlinear dependence (i.e., introducing a variable x n+2n c +1 in place of t in the problem functions and adding the constraint x n+2n c +1 = t) any parametric NLP can be written such that ∇2 xt f (x, t), ∇ t g(x, t) and ∇ t h(x, t) are all constant. 2We consider starting with an approximate solution ( x, λ, μ) which is sufficiently close to the primal-dual solution (x * (t), λ * (t), μ * (t)) at the initial value of the parameter t.Let t be the stepsize for t and consider an iteration of pathfollowing.The predictor-corrector QP is given as, Starting at t = 0 with a given initial approximate solution (x, λ, μ), we iteratively solve (9), take a step in the direction of the primal-dual QP solution, and perform an estimate of the strongly and weakly active sets for t If the original problem satisfies the SSOSC throughout the path, this QP is has a unique solution, as it is strongly convex along the subspace defined by the linearization of the constraints, and this procedure is well defined.
The active sets I + ( x, λ, t) and I 0 ( x, λ, t) are estimated along the path using a technique based on Facchinei et al. (1998), also used in Kungurtsev and Diehl (2014).In particular, for a standard NLP (i.e., n c = 0), we define the optimality residual as and estimate the active sets as follows, where 0 < γ < 1.
For the MPCC pathfollowing algorithms discussed later in Sect.3, we trace the MPCC solutions by tracing some auxiliary parametric NLPs by estimating, at each value of t, the active sets as above, and computing the steps using (9).

Penalty method
In order to obtain a solution of (PMPCC(t)) for a given t, the standard penalty method instead solves the NLP, for some penalty parameter ρ > 0.
There are several important results relating the stationarity properties of the original problem (PMPCC(t)) and its relaxed penalty surrogate problem (P ρ (t)).
Theorem 3.1 (Ralph and Wright 2004, Theorem 5 From the proofs in Ralph and Wright (2004), the minimum ρ that is required for the first two results of this Theorem is equal to, and for the third result, it must satisfy, where ξ is equal to, Note, however, that the MPCC-SOSC does not imply the strong second order sufficiency conditions for the penalty NLP, which is required for the theoretical convergence properties of the predictor-corrector QP (9) for pathfollowing.
One can obtain the multipliers σ * from the solution of the penalty function minimization by taking, where z * are the multipliers for the bound constraints x i ≥ 0, i ∈ N c in the penalty problem (see Coulibaly and Orban 2012).

Description of the penalty pathfollowing algorithm
Now we describe an algorithm that approximately traces the solution trajectory (PMPCC(t)) by following the solution path of the corresponding parametric penalty formulation.
We write the penalty objective as, and correspondingly the Lagrangian of the penalty problem as, with Hessian matrix, We present the Penalty Pathfollowing Algorithm as Algorithm 1.We start with a point close to x * (0) and estimate both the strongly and weakly active inequality constraints using (11).This set includes the complementarity variables treated as constraints, i.e. x i ≥ 0 for i ∈ N c are included and identified as being (strongly) active or not.We shall denote the index set of the complementarity variables deemed strongly active to be N + c (x, (λ, z), μ, t).Note that when discussing the penalty method, η given in (10) includes the dual complementarity measure min(x N c , z), i.e., Note that if the strongly active set is identified correctly and upper level strict complementarity holds, then the solution set of the tightened NLP is equivalent to the MPCC, and thus η is an error bound for the latter.

5:
If Re-estimate the active and the strongly active sets.9: else 10: Set t = t/α.11: end if 12: end while

Theoretical guarantees
The tracking convergence guarantees associated with the penalty algorithm is given below in Theorem 3.2.In particular we use the aforementioned correspondence between necessary and sufficient optimality conditions for the MPCC and the associated penalty problem to show that for a well behaved PMPCC, there is sufficient regularity for the penalty formulation that the solutions are isolated, and thus 1) the problem satisfying complementarity is the only solution to the penalty form and 2) the QP (15), which is its linearization, has a unique solution.Then from standard results in pathfollowing for NLPs the procedure is shown to generate a set of iterates with theoretical tracking guarantees to the parametric solution curve.
In general, for well-behaved problems the procedure successfully traces the solutions, but has several limitations.First, only S-stationary points are guaranteed to be traced.This does not necessarily imply that weaker notions of stationarity that could be spurious, like M or C-stationarity, will not be traced, and it also implies that not all Bstationary points will be traced.Furthermore, an "a posteriori" assumption must hold, that the SSOSC conditions hold for the penalty problem associated with the MPCC.The SOSC for the original MPCC is not enough to ensure the predictor-corrector QP is able to generate the solution path for the MPCC.Theorem 3.2 Assume there exists x * : [0, 1] → R n+2n c such that x * (t) is an isolated S-stationary point for (PMPCC(t)) satisfying MPCC-SOSC and MPCC-LICQ3 for all t, and partial strict complementarity holds for all dual vectors (λ * , μ * , σ * ) corresponding to S-stationarity with x * (t) for all t ∈ D for a set D ⊆ [0, 1] satisfying {0, 1} ∈ D, and [0, 1] \ D is a finite set of points.Then there exists a ρ sufficiently large such that, for all t ∈ D, x * (t) is a solution for (P ρ (t)) at t. Furthermore consider the iterations of Algorithm 1 performed with this ρ appearing in the subproblem (15) with such sufficiently large ρ.Denote by T the set of pathfollowing parameters {t 0 = 0, t 1 , t 2 , . ..}.If it holds that, 1. there is an initial estimate x 0 sufficiently close to x * (0), 2. the SSOSC holds for (P ρ (t)) for x * (t) for all t ∈ [0, 1], 3. there exists only a finite {t a k } ⊆ [0, 1] for which i.e., the number of optimal complementarity active set changes across x * (t) is finite, 4. t k ∈ D for all k, Then it holds that Algorithm 1 successfully generates an approximate solution path {(t k , x k )} with T = {t 0 = 0, t 1 , . . ., t k , . . ., t K = 1} finite, i.e., for any , there exists a 0 such that if Proof Outline: There are two primary steps of the proof.First is to show that there is a one to one locally unique correspondence between x * (t) as a solution of the original problem and as a solution to the penalty formulation.After establishing this, we can use the results on the properties of the penalty formulation above to infer that corresponding strong regularity conditions are satisfied for the penalty problem along the entire path T .What remains is to show that the required minimal ρ for the aforementioned equivalence to hold must be bounded across the homotopy, which corresponds to the quantity defined in (12) having a uniform bound across t ∈ [0, 1].Strong regularity for the penalty problem already implies local boundedness of the primal-dual solution for t ∈ D, thus the onus is to show that the quantity in (12) cannot grow without bound as one approaches some t ∈ [0, 1]\D, which correspond to points at which strict complementarity is lost, i.e., active set changes for the complementarity variables.We split the proof into two cases as depending on the form of the bifurcation occurring at these exceptional points.
Proof First consider x * (0) as given by the conditions of the Theorem, i.e., satisfying S-stationarity for (PMPCC(t)) at t = 0. From Theorem 3.1 we can conclude that x * (0) is a solution to (P ρ (t)) at t = 0 for all ρ ≥ ρ0 for some ρ0 .Furthermore, the LICQ and, by the second assumption in the Theorem, the SSOSC hold at x * (0) as a solution to (P ρ (0)).We can extend this to say that these conditions hold for all x * (t) as solutions to (P ρ (t)), due to the stability of the LICQ and SSOSC (Kojima 1980, Corollary 5.6) and the assumption of x * (t) being an isolated solution path.The final statement of Theorem 3.1, that there is a neighborhood B(t) of x * (t) such that every stationary point of (P ρ (t)) in B(t) is an S-stationary point for (PMPCC(t)), implies there is a one-to-one correspondence between x * (t) as a solution to (P ρ (t)) and for the MPCC in a region around x * (t) for all t ∈ D. Now we show that for all t ∈ D, the right-hand side of ( 12) is uniformly bounded by some ρ for all primal-dual complementarity solutions x * (t), σ * (t), i.e., that the required ρ to satisfy the correspondence between solutions to (PMPCC(t)) and (P ρ (t)) has a uniform lower bound.Now for t / ∈ {t a k } the solution x * (t) is also a solution to the tightened NLP, and since the LICQ and SSOSC hold, its primal-dual solution set is Lipschitz stable (Bonnans and Shapiro 2013, Theorem 4.51).
Together with no components of x * (t) changing from positive to zero (inactive to active) or vice versa for such t, we obtain for any small δ > 0, the uniform boundedness of the quantity in (12) for solutions associated with t ∈ (t a k−1 + δ, t a k − δ).Thus we must only consider the change in ρ1 in (12) across and close to active set changes t a k .Without loss of generality, consider that there is only one complementarity constraint, i.e., n c = 1.Now, given a t c ∈ {t a k }, again without loss of generality, we can consider two exhaustive cases.
Case 1: For t with t c − δ < t < t c , the index i ∈ N c , and, without loss of generality, assume i ∈ N − c , corresponds to an inactive primal solution component (i.e., x * i (t) > 0), and for t satisfying t c ≤ t ≤ t c + δ, x * i (t) = 0 (is active).In addition, x * i+n c (t) = 0 for all t with t c − δ < t < t c + δ.By complementarity of the solution set and Partial Strict Complementarity for t ∈ [t c − δ, t c + δ] \ {t c }, we also have σ * i (t) = 0 and σ * i+n c (t) > 0 for t with t c − δ < t < t c .Given the continuity of the primal-dual solution set to the tightened NLP for t < t c as established in the first part of the proof, it must be that x * i (t c ) = 0 (by the case assumption) and σ * i (t c ) = 0 (by continuity of the tightened NLP multipliers) for a primal-dual solution to the tightened NLP at t c based around x * (t c ).With σ * i (t), σ * i+n c (t) ≥ 0, one can see the one to one correspondence of the dual solution set to the tightened NLP and the set of multipliers satisfying S-stationarity for (PMPCC(t)) at t c with x * (t c ).Now, for t ≥ t c , {i, i and so by the definition of ρ in (12), indices i and i + n c are not included in the set under consideration, and ρ is simply 1 for t c ≤ t ≤ t c + δ.
Thus, for t c − δ < t < t c + δ, we need only be concerned with t < t c , and since σ * i (t) = 0 for such t, we are left with showing boundedness of the quantity and thus the expression i is now unbounded.However, notice that in this scenario, given that x * (t) with σ * (t) as multipliers for the conditions x j ≥ 0, j ∈ N c , is a primal-dual solution to the tightened NLP, we notice that for (P ρ (t)) for any ρ ≥ 0, where the final equalities come from weak stationarity, assuming without loss of generality that there are no general inequality and equality constraints, or that the sum of their multiplier weighted gradients is equal to zero.From this, it can be seen that the necessary optimality conditions of the penalty problem hold for x * (t) with for all ρ ≥ 0. We have already established the local uniqueness of the solutions to (P ρ (t)), thus, this primal-dual pair is the only solution and must be the one the algorithm tracks.
Case 2: it holds that x * i (t) > 0 and x * i+n c (t) = 0 for t < t c , and x * i (t) = 0 and x * i+n c (t) > 0 for t > t c .Note that a multiplier tuple satisfying S-stationarity with x * (t) must satisfy σ * i (t) = 0 for t ≤ t c and σ * i+n c (t) = 0 for t ≥ t c .Thus Partial Strict Complementarity does not hold at t c .Otherwise, for consistency with the assumption on D it must hold that σ * i (t) > 0 for t > t c and σ * i+n c (t) > 0 for t < t c .By the Theorem assumptions, the Algorithm never attempts to solve the pathfollowing problem for t k = t c .At t = t c the situation is now defined exactly as in (Bouza et al. 2008, Proposition 4.4).In particular, for the two tightened NLPs, with x i = 0 fixed for t > t c and x i+n c = 0 fixed for t < t c , x * (t) is optimal with the associated σ * (t) as optimal multipliers for the complementarity constraints now as variable bound inequality constraints.As according to Bouza et al. (2008, Proposition 4.4) there is one homotopy path of solutions that is feasible and it is characterized by, which imply the uniform boundedness of the quantity defined by (12) as t → t c .Now we can bring x 0 as close as necessary to x * (0) such that the conditions of Kungurtsev and Diehl (2014, Theorem 4.2) are satisfied as applied to (P ρ (t)). 4hus the predictor-corrector QP (15) associated with the appropriate active and strongly active set estimates, is able to successfully trace a primal-dual solution with For illustration of the main ideas in the proof, we point to Fig. 1.We now discuss the assumptions in the previous Theorem.The first is standard and natural.The second indicates the gap between the second order necessary and sufficient conditions can only be assumed away, and cannot be guaranteed by any preceding assumptions.The third assumption simply indicates that there are a finite number of Fig. 1 Illustration of the cases described in Theorem 1.The figures on top are potential primal-dual solution paths, as the parameter t changes, the solution undergoes an active set change at t c .This is reflected in the blue lines showing the value of the primal solution of a complementarity pair.The green lines show two possible scenarios of the change in the dual solutions, for Case 1, on the left, and Case 2, on the right, respectively.The figures at the bottom represent the value of thereshold ρ required in the main result, as the parameter t changes.Note that in the first case, this value approaches infinity asymptotically as t → t c , but the structure of the problem remains such that a smaller penalty parameter is adequate active set changes along the solution path, which is expected to hold except for cases of increasingly oscillatory behavior (e.g., with a function of the form sin(1/x)).
It could be of concern that the assumption in the Theorem that t k ∈ D for all k is a posteriori, i.e., depending on the outcome of the Algorithm rather than an a priori condition on the problem itself.However, note that since [0, 1] \ D is finite, also by assumption, and yet t k ∈ R, then this should reasonably be satisfied, and otherwise can be assured with probability one by adding a tiny perturbation to t k .

Active set index method
The next pathfollowing algorithm is intended to trace B-stationary points, which as indicated describe the tightest computable notion of necessary first order optimality.Their importance and the combinatorial approach of the active set algorithm described here is inspired by the "SLPEC-EQP" algorithm presented in Leyffer and Munson (2007) and more recently (Kirches et al. 2022) [see also Lenders (2018, Chapter 7)].
In this line of work, it was observed that M, C and W-stationary points could have linearly feasible descent directions, and B-stationarity is the precise stationarity concept associated with the lack of such descent directions.S-stationarity is the strongest notion of stationarity, but an algorithm with guarantees in finding only S-stationarity could neglect finding points that are B-but not S-stationary.
Consider the LPCC, where β > 0 is a trust-region parameter.It is clear from considering the LPCC associated with the definition of B-stationarity ( 6) that d = 0 being the solution to ( 16) corresponds to B-stationarity.The algorithms in Leyffer and Munson (2007) and Kirches et al. (2022) solve this LPCC, using the solution to estimate the active set and subsequently perform acceleration steps on the active set.The authors prove the global convergence of this procedure to a B-stationary point.Recall that d = 0 solving the LPCC (( 6) and ( 16)) at some (x * , t) implies that every index I NLP for I ∈ J (x * , t) is optimal.We focus on estimating the set J (x * , t), noting that the double active set indicates possible branching in the directional derivative of the solution.
Thus, given t and an initial point x that we expect to be close to a solution for problem (PMPCC(t)) at t, we first seek to estimate the active set of complementarity variables at the optimal point closest to x. Denoting B(t) as the set of B-stationary points for (PMPCC(t)) at t, we write this as, Next we must identify the set of index I NLPs which requires a reliable active set identification procedure for MPCCs.However, there is no certifiably accurate method to do this and we make the following assumption.
Assumption 1 We have a reliable means of estimating A 0 (x, t) from x, if x is sufficiently close to x * (t), a B-stationary point of (PMPCC(t)).
Heuristically, we consider i as active, and thus a potential member of the index I NLP when x i ≤ 0 for some (problem-tuned) tolerance 0 .Practically, we form the set of branches based on the estimate of the corresponding doubly-active set, i.e., {i : We then use this to form an estimate of the set J (x * , t) and perform parametric pathfollowing on every NLP.Denote this estimate to be C 0 (x, t), i.e.,

Description of the algorithm
We present the Active Set MPCC pathfollowing algorithm as Algorithm 2. First, given an approximate B-stationary solution x to the MPCC at t = 0, we estimate and prune the set of index I NLPs.For each I ∈ C 0 (x, 0) defined in ( 17), we form the NLP and solve it using an appropriate solver for standalone NLPs in order to obtain the multiplier estimates for the initial (PMPCC(t)) at t = 0.If no appropriate dual variable approximate solutions (i.e., with the right signs) exist, then this I is not in C 0 (x, t) and we discard it.Otherwise, for every I in C 0 (x, t), the corresponding predictorcorrector QP is then used to obtain the next solution.The QP we solve is the following (assuming, again, affine as welldependence of g(x, •) and ∇ x f (x, •) with respect to the parameter t).
where A 0 (I ) is the active set designated by index I .We solve the QP for every feasible I to obtain a set of solutions { x I }.
We then check, after the step, if the solution still satisfies the full set of complementarity constraints.Then, for each I , we recalculate C 0 (x + x I , t + t) if necessary and investigate if the new primal solution estimate remains stationary for every single index I NLP in C 0 (x + x I , t + t) by obtaining the appropriate dual solution minimizing the stationarity norm, min where I + ⊆ N c is a set such that i ∈ I + implies x i = 0.In particular, if the objective value of ( 19) is less than 0 then we consider it approximately stationary in our Algorithm.
If the point x + x I is not feasible, or not approximately optimal for every I ∈ C 0 (x + x, t + t), then we cut this branch and proceed to starting the next branch, i.e., move on to the next I .Finally, if it tracks the solution until the final iterate at t = 1, for the sake of completeness we pathfollow the other branches as well.

15:
If there is more than one such x I , in parallel branch the algorithm from each new primal solution estimate.

16:
If for all I , x I = 0, set t = t/α.17: Re-estimate the double complementary active set C 0 (x + x I , t + t) by ( 17).18: end while

Theoretical guarantees
Given that the algorithm consists of pathfollowing the set of index I NLPs, the convergence properties can be inferred from the corresponding results for standard NLP pathfollowing using the predictor-corrector QP as presented in Kungurtsev and Diehl (2014).
Assume this set is maximal, i.e., there is no K > K that strictly includes {x * k } k=1,...,K such that these conditions are satisfied.Given an initial point x 1 (0) sufficiently close to x * k (0) for some k, and desired threshold D, Algorithm 2 will generate a set of iterates {x l (t j )} l=1,..., K , j=1,...,J l with K ≤ K such that x l (t j ) − x * k (t j ) ≤ D for some k ∈ {1, . . ., K } for all l ∈ {1, . . ., K } and j ∈ {1, . . ., J l } and t j ∈ [t a k , t b k ] for this k and max j {t j } = 1.
Proof Under Assumption 1, since B-stationarity is equivalent to a set of NLPs satisfying the first order necessary optimality conditions, x 1 (0) is an approximate stationary point to each Index I NLP found in C 0 (x, 0).If there exists a branch x * k (t 0 + t) for any t > 0, then by Kungurtsev and Diehl (2014, Theorem 4.2) the active set estimation and predictor-corrector QP will trace it to desired accuracy.By the fact that the distance to the solution of any NLP is O(•) of the value of the optimality residual η(•) (see, e.g.Hager and Gowda 1999), a true B-stationary point will satisfy a low objective value of (19) for every index I NLP.
If during the pathfollowing another constraint becomes doubly-active at t, it is added to C 0 (x(t), t), and subsequently traced from t on.

Examples
In this section we present the results of Algorithms 1 and 2 on a set of toy problems and a phase transition engineering case study.We performed the experiments using MATLAB 2018a.The initial point was obtained using fmincon with the SQP option, constraint tolerance 1e-7 and step tolerance 1e-12 for the penalty algorithm, and the multistart based global optimization toolbox for the active set algorithm.Quadprog with the default settings was used to solve the subproblems.A constraint is considered active if it is so approximately under a tolerance of 1e-5.t is initialized to 0.1 and is increased or decreased by a factor of α = 1.5 upon a successful or unsuccessful pathfollowing step.

Example 1
Consider the problem, For t < 0, (0, −t) is a local minimizer and the origin is W-stationary.For t = 0, the origin is a strongly stationary local minimizer.For t > 0, (t, 0) is a local minimizer and the origin is W-stationary.

Penalty Method
The penalty method is able to successfully trace the solution (0, −t) The primal solutions as well as the complementarity multiplier estimates given by ( 13) are shown in Fig. 2.

Active Set Method
The active set method is able to successfully trace the solution (0, −t) for t ∈ [−1, 0] through (t, 0) for t ∈ [0, 1].It recognizes the double-active set at t = 0 and discards the branch of the Index I NLP associated with the path of solutions for t ∈ [−1, 0], switching to the other branch thereafter.

Example 2
This problem is inspired by Leyffer and Munson (2007).
Fig. 2 Illustration of the traced solutions to Example 1 found using the penalty algorithm.The circles are the points visited by the method.As expected, since the solutions are all S-stationary, the complementarity multipliers are zero when the corresponding variable is inactive and become positive as they become active, and vice versa.For t < 0, the origin is the local minimizer and only stationary point.For t > 0, the point (t, 0) is a local minimizer and the origin is M-stationary.

Penalty Method
In this case, the penalty method always traces the origin for t ∈ [−1, 0], and the local minimizer (t, 0) for t ∈ [0, 1].Thus it does not appear to follow the M-stationary point.The set of points traced along the homotopy is shown in Fig. 3.

Active Set Method
The active set method traces two solutions along the bi-active set for t ∈ [−1, 0].Subsequently, the solution branch for which x 1 = 0 stops pathfollowing at t = 0 and the one with x 2 = 0 pathfollows successfully along (t, 0) for t ∈ [0, 1].The solutions are shown in Fig. 4.

Example 3
This problem arises from Scheel and Scholtes (2000).Consider, Fig. 5 Illustration of Example 3. The level sets of the objective function at t = −1 and t = 1 are drawn as circles.As t traces from −1 to 1, the level sets move continuously from the bottom left to the top right.The local minimizers are marked as red dots.Up until t = 0, i.e., for t ∈ [−1, 0], the unique local minimizer is the origin.For t ∈ (0, 1] there are two local minimizers, (t, 0) and (0, t).In addition, there is a C-stationary point at the origin, for which clearly both axes are a descent direction, indicating the weakness of this stationarity condition For t ≤ 0, the origin is the unique minimizer, and is S-stationary.For t > 0, the origin is C-stationary, and (t, 0) and (0, t) are both local minimizers.For any t, the Lagrangian Hessian is ∇ 2 x x L = 2 0 0 2 , as this matrix is positive definite, it is also positive definite on any cone of directions.We illustrate the problem in Fig. 5.

Penalty Method
Note that for t = 0, the Lagrangian Hessian for (P . At the origin, every multiplier is null, i.e., σ * 1,2 = z * 1,2 − ρx * = 0, thus for (P ρ (t)) at t = 0, no constraint is strongly active and both variables are weakly active at the origin, i.e., Thus, although ∇ 2 x x L p is positive definite along the cone of feasible directions, which are 1 0 T and 0 1 T , for proper pathfollowing of (P ρ (t)) using ( 15), the strong SOSC conditions must hold for (P ρ (t)) at x * , and ∇ 2 x x L p must be positive definite on the subspace orthogonal to the null-space of the strongly active constraint gradients, which is trivially R 2 , and thus includes 1 −1 T .Thus, predictably the penalty algorithm is not able to pathfollow past t = 0, in particular, as the reduced Hessian in the constraint gradient null-space of the pathfollowing QP ( 9) is indefinite, the corresponding QP is nonconvex.Thus, theoretical guarantees on NLP pathfollowing do not apply to the corresponding predictor-corrector QP for this problem.
The optimality conditions for the QP (9) based from the origin are, From these equations, it must hold that x 1 = x 2 > 0. Note that this is not a feasible direction for the original problem, and thus the QP for the pathfollowing of (P ρ (t)) for any ρ > 0 does not correspond to pathfollowing for the original MPCC.The penalty method is able to trace the origin solution for t < 0, and for t > 0 it appears to pick either local minimizer with equal probability, while never converging to the C-stationary point, and finds it impossible to pathfollow too close to t = 0.

Active Set Method
The active set method, started at t = −1, recognizes the presence of the doubly-active set and proceeds with two separate solution paths, one with fixed x 1 = 0 and the other with fixed x 2 = 0.The algorithm is able to follow the origin for t ∈ [−1, 0] and each path follows the separate local minimizers ((t, 0) and (0, t) for t ∈ [0, 1]).The two branches of solutions are shown in Fig. 6.

Example 4
This problem is original to this paper.It illustrates how the penalty method can fail to find local minimizers that are B-stationary but not S-stationary, whereas the active set method can do so.Note that this problem fails the MPCC-LICQ condition.
For all t ∈ [0, 1] the origin is the unique local minimizer.The stationarity conditions for the problem are, ) is M-stationary, and furthermore for the index {1} NLP, the second set of multipliers is optimal as well as (0, −2, 2t), and for the index {2} NLP the first set is optimal, and thus the origin is also B-stationary for all t ∈ [0, 1].

Penalty Method
The standalone penalty method is unable to find a satisfactory initial point.Initializing from one of the stationary points does not result in any successful pathfollowing step.

Active Set Method
The active set method is able to trace both index I NLPs for this set of problems, tracing the origin for t = [0, 1] along the multipliers (0, −2, 2t) for the index {2} NLP and (2, 0, 2t − 2) for the index {1} NLP.The primal and dual solutions traced are shown in Fig. 7.

Example 5
This problem arises in the MacMPEC collection of MPCC problems (Leyffer 2000).

Penalty Method
The penalty method is unable to find a satisfactory initial point, i.e., a solution at t = 0. Initializing from one of the stationary points does not result in any successful pathfollowing step, with the error that every predictor-corrector QP problem is nonconvex (note that the problem itself is linear, but the penalty method adds a bilinear term to the objective).For t ∈ [ 3 4 , 1] the algorithm is able to pathfollow the solution x * = (0, 0, 0).

Active Set Method
The active set method pathfollows along the origin, doing so twice, once for each index in the double active set, for t ∈ [0, 1].

Example 6
We now consider a nonlinear example, to demonstrate the properties of the pathfollowing algorithms on a nonlinear problem.
The constraints and solutions are shown in Fig. 8a.The point must lie outside of the circle and to the left of x 1 ≤ 1, and the objective function seeks to increase x 1 and decrease x 2 .As t changes from −1 to 1, the circle expands, modifying the feasible region to exclude any part of the x 1 -axis.The solution is always S-stationary.We plot the result of the penalty algorithm applied to this problem in Fig. 8b and the active set based procedure performs similarly.

Flash calculation example
This section describes a classic chemical engineering problem appearing in, e.g., Skogestad (2008, Chapter 7).It models vapor/liquid equilibrium, wherein the Gibbs energy is minimized at some given temperature and pressure.In this case we consider pressure to be fixed and vary the temperature.Pathfollowing for a similar PMPCC was considered in Nakama et al. (2022), using a different approach that does not attempt to trace each branch.
During the operation of chemical processes, sometimes the thermodynamic phase in a system is not known a-priori.It may be liquid, gas, or a mixture of both, and each of these three cases must be described by a different set of equations.For example in natural gas liquification plants, some sections of the heat exchangers may be filled with liquid, which evaporates as the operating conditions change.A priori, it is not known what conditions will be present, and when optimizing the operation one can model the system using complementarity constraints.A large scale optimization model of a full process is out of the scope of this article, but we present a case as a proof of concept that the algorithms from this paper can be applied successfully to such a system.
We consider an exogenous feed stream F that is separated into vapor V and liquid L product, of three different compounds {y i } and {x i }, respectively.The relative split of F to V and L as well as the quantity of each of the three compounds in each phase are considered variables we must solve for.They must satisfy i x i = 1 and i y i = 1, and x i L + y i V = z i F for the total given quantities {z i }.Using Raoult's law, we have that, where p sat i is a nonlinear function of temperature, where A i , B i and C i are physical constants. 123 The fraction a ∈ [0, 1] = V /F of gas is determined by the Rachford-Rice equation, however, with a caveat.Specifically, attempt to solve the following equation for a t : Now, the bubble point temperature is T b = 382.64 and the dew point T d = 393.30,and thus for T < T b , the solution for a t will be negative, but clearly the ratio a cannot be negative, instead it implies that every substance must be a liquid, i.e., a = 0 and similarly for T > T b , the solution for a t will be greater than one, but clearly every substance must be gas, i.e., a = 1.Thus, we need to, after solving the Rachford-Rice equation for a t , obtain the value of a by max(0, min(1, a t )), a nonsmooth operation problematic for gradient-based optimization solvers.To avoid this, the true ratio of gas a will be determined by the complementarity system, where we introduce slack variables s l and s v complementary with L and V (i.e., s l is zero if L is positive, and L is zero if s l is positive) to compute a from a t by solving the first equation.
The resulting problem is a nonlinear complementarity system involving the 13 variables V , L, {x i } i=1,2,3 , {y i } i=1,2,3 , a, a t , T , s v and s l .In addition, the nonlinear equations relating temperature to the component distributions are ill-conditioned.Thus, to define the MPCC, we, • lift Raoult's law and the Rachford-Rice equation, making K i , 1 (K i −1) (which we write as k i ), log p sat i variables, alongside V , L, x i , y i , a, a t , T , s v and s l for a total of n = 22.
• We consider T as both a parameter and a variable, by defining an additional equality constraint T = T target , and we trace the solution for T target = 380 to 400.• If we were to add the constraint defining the proportial of gas as a = V F , then all of these equations define a nonlinear complementarity problem (NCP).However, when taking into account active constraints at the solutions of the NCP, the entire system does not satisfy MPCC-LICQ, thus we do not make the vapor fraction a = V /F a hard constraint, but define our objective to be, and thus define an M PCC.In summary, the MPCC we solve is defined as follows, min V ,L,{x i },{y i },a,a t ,T ,s v ,s l 1 2 (a F − V ) 2 s.t.k i = 1 (K i −1) , i = 1, 2, 3 q i = log p sat i , i = 1, 2, 3 q i = log(10)(A i − B i T +C i ), i = 1, 2, 3 The constants we use in the experiments are given in Table 1.

Penalty method
The penalty method accurately tracks the solution, achieving the liquid-gas balance as expected.
Every solution is a S-stationary point, and the objective function is strictly convex, which means that for any point satisfying the first order necessary conditions for optimality, the second order sufficient conditions also hold.Thus all the necessary conditions for the convergence of the penalty method for solving the MPCC and tracing a parametric curve are satisfied.
The solutions for all of the variables are given in Fig. 9.As expected, the variables X and s x , as well as Y and s y are complementary, with the slacks being positive when their corresponding complementary variables are zero, and zero when they are inactive.The proportion of gas a corresponds to V .The variables k i (representing 1 (K i −1) ) appear to vary the most nonlinearly, resembling an exponential function with respect to the parameter, which is equivalent to the temperature.
Note that the required step-size for pathfollowing is generally fairly small, and gets smaller during active set changes.

Active set method
The active set method also traces the set of solutions.At each active set change, a double-active set was recognized.The previous branch was eventually cut, and the other branch initiated.Although the method converged for this problem, it took considerable tuning with regards to tolerances for the double-active constraints, on the one hand, and the required tolerance for checking the other index I NLPs, while tracing a branch, on the other.In particular, numerically, a double-active constraint will stay active for a few homotopy steps before being seen as no longer bi-active, and meanwhile the steps must still be accepted as valid for all index I NLPs.

Discussion
The two phase flow example illustrates the general efficacy of the two Algorithms on a real problem.The simpler problems indicate some subtle distinctions of the two methods.In particular, the penalty method is easier to tune and implement, traces S-stationary but not B-stationary points, and also can trace C-stationary points, which may not be desired by the user.By contrast, the active set method is able to trace B-stationary points and does not trace spurious stationary points.However, since it must trace every combination of selections of doubly-active constraints, it is generally slower and more difficult to tune for problems with many complementarity constraints and suspected bi-activity.

Conclusion
In this paper we studied the parametric mathematical program with complementarity constraints, presenting the first investigation of algorithmically pathfollowing these programs.We formulated two algorithms, one based on the penalty method for standalone MPCCs, and one based on following the possible branches of NLPs suggested by double-active constraints and the B-stationarity condition.We studied their performance on a set of examples constructed to distinguish the types of parametric properties there could be, illustrating the performance on this set of problems, as well as a "real-world" problem associated with chemical process engineering.
Pathfollowing parametric optimization problems can be important in a number of applications, including nonlinear model predictive control, as well as assessing the practical sensitivity of a solution with respect to a parameter in industrial design.We intend that this work presents a solid first step in the development of literature and algorithms for the pathfollowing of parametric MPCCs.In general, the tradeoffs associated with seeking tighter versus looser notions of stationarity inherent to standalone methods for solving MPCCs carries over to pathfollowing parametric ones, and we expect this to enrich further understanding of MPCCs in general and algorithms for pathfollowing parametric MPCCs in particular.

Fig. 3
Fig. 3 Illustration of the penalty algorithm solutions for Example 2

Fig. 4
Fig. 4 Illustration of the active set algorithm solutions for Example 2

Fig. 6
Fig. 6 Illustration of the active set algorithm solutions for Example 3

Fig. 7
Fig. 7 Illustration of the active set algorithm solutions for Example 5 Fig. 8 Nonlinear Example

Fig. 9
Fig. 9 Primal solutions for the Liquid-Phase flow problem .1-5.2, Proposition 5.3) • If x * is S-stationary for (PMPCC(t)) at t then for ρ sufficiently large, x * is stationary for (P ρ (t)).Furthermore, if the MPCC-LICQ, MPCC-MFCQ, or MPCC-SOSC hold for (PMPCC(t)) at t at x * then the LICQ, MFCQ or SOSC hold at (P ρ (t)) at t at x * , respectively.• If x * is stationary for (P ρ (t)) and x * ∈ X (t), then it is S-stationary for (PMPCC(t)).Furthermore if the LICQ and SOSC hold at x * for (P ρ (t)) at t, then the MPCC-LICQ and MPCC-SOSC hold at x * for (PMPCC(t)) at t. • If x * is S-stationary with associated (λ * , μ * , σ * ) for (PMPCC(t)) at t, and the MPCC-LICQ and partial strict complementarity hold at x * , then for ρ sufficiently large, there is a neighborhood B of x * such that every stationary point of (P ρ (t)) in B is an S-stationary point for (PMPCC(t)).

Table 1
Constants used in the flash calculation