# ParametricNDSolve KKT Constraint Directly

I am trying to directly solve for an optimal control using Pontryagin Maximum, where H is the Hamiltonian, in the presence of a KKT constraint. The constraint is ignored. I am wondering if a direct solution is impossible due to an inability to implement the necessary KKT requirement that lambda3 >= 0. The principle is whether ParametricNDSolveValue, or a surrogate, be made to respect the inequality constraint? The result of this question will be of considerable interest to economists, chemical engineers and others.

uu = E^(-0.015 t) (1 - c[t]^-1);
constraint1 = k'[t] == inv[t] - 0.04 k[t];
constraint2 = inv[t] == k[t]^0.85 - c[t];
inequality = 5 10^6 - k[t] >= 0;
H =
uu + lambda1[t] constraint1[] +
lambda2[t] (constraint2 /. Equal -> Subtract) + lambda3[t] inequality[];
eqn1 = D[H, k[t]] == -lambda1'[t];
eqn2 = D[H, c[t]] == 0;
eqn3 = D[H, inv[t]] == 0;
model = {eqn1, eqn2, eqn3, constraint1, constraint2, lambda3[t] inequality[] == 0};
statevars = Union[Cases[model, _Symbol[t], Infinity]];
trial = {1.8, 7};
testinitialvals =
Solve[model /. t -> 0 /. Thread[{c, k} -> trial]][]
With[{pv1 =
ParametricNDSolveValue[{model, c == c0, k == k0}, statevars,
{t, 0, 100}, {c0, k0}, MaxSteps -> 10^8]},
solver[params_] := Thread[statevars -> Apply[pv1, params]]]
solver[trial];
{Plot[k[t] /. %, {t, 0, 100}], Plot[inequality[] /. %, {t, 0, 100}]}


@Daniel - customising equations might achieve a result in specific cases such as this toy problem. Mathematica 10 Documentation (ref/DiscreteVariables Applications) has some useful examples and there is a good example of equation switching here.. In this toy problem most of this switching leads to an undefined solution in some of the variables. I am interested in a general purpose approach for DAE systems with hundreds of equations and constraints where non-systematic hand-crafted customization is not feasible. It would have been useful if something like WhenEvent[lambda3[t] < 0, {a[t] -> lambda3[t], lambda3[t] -> 0}] (i.e. lambda3 goes to 0 from above) was way of implementing the lambda3 >= 0 constraint. While I am no expert with WhenEvent, it seems to me that this doesn't not work unless there is some other factor that will cause ParametricNDSolveValue to pick-up on lambda3's transition from zero to non-zero when the associated constraint becomes active.

@Jens & OP - thanks for placing my response to Daniel back in the question

Update 1: I find this problem quite difficult to solve with the discounted objective function and more than just one equality and one inequality constraint. Also, implementing a WhenEvent approach is difficult because the right-hand-side of the model equations cannot be switched (since there are different equations in each case) and NDSolve complains of being underdetermined with a lesser set of equations.

However, here is a formulation for a slightly simplified version of the problem that does work. It demonstrates the correct action of the inequality constraint using the two solved model solutions:

uu = Limit[(c[t]^(1 - \[Theta]) - 1)/(1 - \[Theta]), \[Theta] -> 2];
constraint = k'[t] == k[t]^0.85 - c[t] - 0.04 k[t];
kmax = 5 10^6; inequality = kmax - k[t] >= 0;
H = uu + lambda1[t] constraint[] + lambda3[t] inequality[];
eqn1 = D[H, k[t]] == -lambda1'[t];
eqn2 = D[H, c[t]] == 0;
(*Solution by discretizing the inequality's prebinding & postbinding*)
model = Solve[{eqn1, eqn2, constraint,
lambda3[t] inequality[] == 0}] /. Rule -> Equal;
Union[Cases[model, _Symbol[t], \[Infinity]]];
statevars = Join[%, D[%, t]]; trial = {1.8, 8}; testt = 100;
prebinding =
statevars, {t, 0, testt}];
tbinding = t /. FindRoot[k[t] == kmax /. prebinding, {t, testt}];
bindingvals = Flatten[{c[t], k[t]} /. prebinding] /. t -> tbinding;
postbinding =
NDSolve[{model[],
Thread[({c[t], k[t]} /. t -> tbinding) == bindingvals]},
statevars, {t, tbinding, testt}] // Quiet;
result = Map[# ->
Piecewise[{{# /. prebinding, t < tbinding}}, # /. postbinding] &,
statevars];
inv[t] = k[t]^0.85 - c[t];
Map[Plot[# /. result, {t, 0, testt}, PlotLabel -> #] &,
Join[statevars, {inv[t]}]]


The "pure" inequality constraint in this toy problem is fairly simple. It seems to me that this method is not scalable for non-toy problems with many "mixed" inequality constraints (i.e. comprising both state and control variables in nonlinear equations) because (a) the much more complicated models will not be amenable to neat solution; and (b) the large number of mixed constraints (say 70) may bind and not-bind at various times creating a huge number of potential sets of equations.

Perhaps NDSolve could easily handle this task if it was possible to specify some output variables as >=0 (i.e. a critical issue for the KKT lambda multipliers). This is indeed a possibility in the underling Sundials IDA solver but not currently available through the NDSolve parser. One possibility is using the Mathematica interface for R and the R interface for Sundials. Another possibility might be to use an SQP solver.

Update 2: Above I suggest that another possibility might be to use an SQP solver. This is a different paradigm than that proposed in my original question. An SQP solver is one option for a "discretize then optimize" procedure. Many people are not aware that "any old NLP solver" is not suitable for this type of problem. Bi-level programming is another type of algorithm used. Both SQP and bi-level have had some success for local solutions, albeit limited and hard-won according to the literature. Also, this procedure doesn't scale well with say 200 variables, 1,200 time steps and Runge-Kutta 7th order accuracy!

The heart of my original question (above) is the reverse procedure, viz "optimize then discretize." This is closest to optimal controls theory and would achieve a clear symbolic optimum for rapid numerical solution by NDSolve (i.e. discrete Runge-Kutta).

My Update 1 showed some tantalizing prospect that an "optimize then discretize" approach with NDSolve may be possible. There has been minor additional success using ImplicitRegion to set a lambda>=0 bound. However, it doesn't work with the k[t] inequality, although it does with an an additional constraint i.e "350000 - inv[t] >=0," as long as "lambda4[t] x (350000 - inv[t]) = 10^-6" instead of 0 (which is one work-around used for Mathematical Programs with Equilibrium Constraints (MPEC) optimal controls).

The primary issue is that MPECs are very hard to solve. This is because KKT inequality constraints with the complementary condition "lambda x inequality = 0 && lambda >= 0 && inequality >=0" do not prima facie satisfy KKT's own constraint independence requirements (Linear Independence LIQC & weaker Mangasarian-Fromovitz MFQC Constraint Qualifications)! After 2 decades of research a few MPEC solvers for complementary constraints are only now emerging. However, the techniques have yet to make their way into industrial-grade DAE solvers such as NDSolve.

While I need MPEC capability in NDSolve for elegant rapid optimal controls solutions, I have developed a snazzy "discretize then optimize" approach that doesn't suffer too many of the disadvantages of traditional solvers like SQP. My approach uses Mathematica functions, defined by handles, for the optimal controls. NMinimize calls ParametricNDSolveValue to solve the DAE at each iteration. ParametricNDSolveValue is incredibly fast as only the parameters are in the iteration loop. The "discretize then optimize" approach has a limited number of optimization variables (just 80 or so handles for a very big problem). Also, inequality constraints are a breeze. From a technical standpoint the objective function doesn't need to be "twice-differentiable" as analytically required for the "optimize then discretize" approach. So there is much more flexibility in objective functions for multi-goal programming. I am writing a paper on this method, which I will add here when ready.

The answer to my original question does appear to be "impossible" for now but I hope perhaps this situation will be remedied within a year or so.

• I'm not convinced about this closing. The question is very much connected to Mathematica computations. They are difficult, granted. A basic restatement is this: Can ParametricNDSolveValue, or a surrogate, be made to respect the constraint? I don't know the answer but I think the question at least fits the general scope of this site. Too localized? Maybe. I'm not sure though. Sep 29, 2014 at 15:34

Not a complete answer but might give some ideas. As noted, it is impossible to enforce the inequality since it can conflict with the equations. One possibility might be to cap k[t] and force the derivative to vanish when it hits the cap. Below is code that could be used.

constraint1 =
k'[t] == (inv[t] - 0.04 k[t])*
Piecewise[{{1, 5 10^6 - k[t] > 0}}, 0];


Now we can get rid of the inequality constraint and lambda3.

Another possibility would be to have an event handler mechanism that bails when k[t] gets too large. Which route to take will depend on what you want k to do once it hits that point.

• OP had a response to you in a deleted answer (which is now an addendum to the question). Please respond under the question to notify them (if you choose to respond).
– rm -rf
Oct 2, 2014 at 6:25