5
$\begingroup$

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?

$\endgroup$

1 Answer 1

3
$\begingroup$

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}]

enter image description here

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}]

enter image description here

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.