# Optimal control via ParametricNDSolve

I'm wondering about the possibility of employing ParametricNDSolve to solve a class of constrained optimal control problems.

Here's an example: The system under consideration is $$\begin{pmatrix} \dot x_1(t)\\ \dot x_2(t)\end{pmatrix} = \begin{pmatrix} x_2(t)\\ u(t) \end{pmatrix},$$ with initial condition $$x_1(0) = -1/2$$ and $$x_2(0) = -1$$. I have a specific form of the control trajectory given by $$t\mapsto u(t) = c_1 \frac{1}{\sqrt{\pi\mathcal D}} e^{-t^2/\mathcal D} + c_2 \frac{1}{\sqrt{\pi\mathcal D}} e^{-(t-1)^2/\mathcal D}$$, where $$c_1, c_2$$ are real-valued coefficients. The objective is to minimize $$\int_0^2 u(t)^2 \,\mathrm{d}t$$ subject to the preceding controlled differential equation while ensuring that, for instance, $$x_1(t) \ge -1$$ for all time $$t\in[0, 2]$$.

I begin with the initialization:

tInit = 0; tFin = 2; paramD = 2;

Then I build the solution:

solution = ParametricNDSolve[
{x1'[t] == x2[t],
x2'[t] ==
c1* E^(-t^2/paramD)/Sqrt[\[Pi] * paramD] +
c2*E^(-(t - 1)^2/paramD)/Sqrt[\[Pi] * paramD],
x0'[t] == (c1* E^(-t^2/paramD)/Sqrt[\[Pi] * paramD] +
c2*E^(-(t - 1)^2/paramD)/Sqrt[\[Pi] * paramD])^2,
(* initial conditions *)
x1[tInit] == -1/2, x2[tInit] == -1,
x0[tInit] == 0},
(* declare the states of the ode *)
{x1, x2, x0},
(* declare the time interval *)
{t, tInit, tFin},
(* declare the parameters *)
{c1, c2}];


Here the variable $$x_0$$ keeps track of the objective function; the objective is to minimize $$x_2(2)$$, I think it should be possible to pass it on to NMinimize. For instance, extracting an optimal pair $$(c_1^*, c_2^*)$$ with linear constraints on $$c_1, c_2$$ can be done via:

NMinimize[{x0[c1, c2][tFin] /. solution,
1 <= c1 <= 4 && -5 <= c2 <= -1}, {c1, c2}]


However, this optimal pair $$(c_1^*, c_2^*)$$ solves the minimization problem without the pathwise constraint $$x_1(t) \ge -1\;\forall t\in[0, 2]$$ and with unnecessary constraints on $$(c_1, c_2)$$. Is there a mechanism to include the constraints on $$x_1$$ as mentioned above?

In principle, this problem can be solved by defining the function,

f[c1_?NumericQ, c2_?NumericQ] :=
NMinValue[{x1[c1, c2][t] /. solution, tInit < t < tFin}, t]


And then minimizing

NMinimize[{x0[c1, c2][tFin] /. solution, f[c1, c2] >= -1,
3.5 <= c1 <= 4 && -1.5 <= c2 <= -1}, {c1, c2}]

(* {0.892063, {c1 -> 3.91284, c2 -> -1.29659}} *)


and the corresponding variables are

Through[{x1, x2, x0}[c1, c2]] /. Last[%] /. solution;
Through[%[t]];
Plot[%, {t, tInit, tFin}, AxesLabel -> {t, "x1,x2,x0"},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}]


This computation requires, about 15 minutes. A faster but approximate approach is to plot f[c1, c2] and x0[c1, c2][tFin] /. solution.

Show[
ContourPlot[x0[c1, c2][tFin] /. solution, {c1, 3.5, 4}, {c2, -1.5, -1},
FrameLabel -> {c1, c2}],
ContourPlot[f[c1, c2] == -1, {c1, 3.5, 4}, {c2, -1.5, -1},
ContourStyle -> {Black, Thick}],
ImageSize -> Large, LabelStyle -> {15, Bold, Black}]


We see that the f = -1 line falls furthest below the contour 0.9 neac {3.9, -1.3}. A good approximation can be obtained by evaluating x0[c1, c2][tFin] /. solution along the f = -1 line and then finding its minimum value there.

f0 = ContourPlot[f[c1, c2] == -1, {c1, 3.5, 4}, {c2, -1.5, -1}];
pts = Cases[Normal[f0], Line[z_] -> z, Infinity] // First;
val = Through[(x0 @@@ pts /. solution)[tFin]];
Ordering[val, 1] // First;
{val[[%]], pts[[%]]}

(* {0.892076, {3.91095, -1.29464}} *)


Addendum: Largely Symbolic Solution

The ODEs in the question also can be solved symbolically. For x1[t],

s = Collect[DSolveValue[{x1'[t] == x2[t],
x2'[t] == c1*E^(-t^2/paramD)/Sqrt[π*paramD] +
c2*E^(-(t - 1)^2/paramD)/Sqrt[π*paramD], x1[tInit] == -1/2,
x2[tInit] == -1}, x1[t], {t, tInit, tFin}], {c1, c2}, Simplify]

(* -(1/2) - t + c2 (E^(-(1/2) (-1 + t)^2)/Sqrt[2 π] - 1/Sqrt[2 E π]
- 1/2 Erf[1/Sqrt[2]] + 1/2 t Erf[1/Sqrt[2]]
+ 1/2 (-1 + t) Erf[(-1 + t)/Sqrt[2]])
+ c1 ((-1 + E^(-(t^2/2)))/Sqrt[2 π] + 1/2 t Erf[t/Sqrt[2]]) *)


The corresponding function g, similar to f defined above, is

g[c1t_?NumericQ, c2t_?NumericQ] :=
NMinValue[{s /. {c1 -> c1t, c2 -> c2t}, tInit < t < tFin}, t]


and the objective function is

i = Integrate[(c1*E^(-t^2/paramD)/Sqrt[π*paramD] +
c2*E^(-(t - 1)^2/paramD)/Sqrt[π*paramD])^2, {t, tInit, tFin}]

(* (2 c2^2 Erf[1] + (2 c1 c2 (Erf[1/2] + Erf[3/2]))/E^(1/4)
+ c1^2 Erf[2])/(4 Sqrt[π]) *)


From these quantities, optimal values are obtained from

NMinimize[{i, g[c1, c2] >= -1, 3.5 <= c1 <= 4, -1.5 <= c2 <= -1}, {c1, c2}]

(* {0.892063, {c1 -> 3.91284, c2 -> -1.2966}} *)


in about 15 minutes. If decimal instead of exact expressions are used in the computation,

{sN, iN} = {s, i} // N


the same results are obtained in 8 minutes.

• Thank you; your answer has indicated some rather interesting directions to me. Commented Jul 15, 2022 at 7:17
• My sincere thanks once again for the addendum. Commented Jul 16, 2022 at 3:56