2
$\begingroup$

I have a system of ODEs say

$$x'(t)=2 x(t)-x(t)y(t)-3 x(t) g(t),$$ $$y'(t)= y(t)-2 x(t)y(t)- x(t) g(t)$$

with ICs $x(0)= 1, y(0)=1/2$, say, for $t_0 \leqslant t \leqslant t_f$, where $t_0$ and $t_f$ are start and final times.

The function $g(t)$ is a discrete function: $g(t) = u_i$ for $t_i<= t <= t_{i+1}$ with $t_i$ a discretization of the interval $[t_0,t_f]$, $t_i= t_0+i h$ with $h=(t_f-t_0)/N$, a fixed step size. Note that the parameters $u_i$ are either $0$ or $1$.

The aim is to find the optimal sequence ${u_0, u_1,...,u_{N-1}}$ with $u_i = 0$ or $1$ such that an objective function, say

$$\int_{t_0}^{t_f}[ (x(t)-1)^2+(y(t)-1)^2]dt$$

is minimum.

Any help is very much appreciated. Can try with $N=5$.

Many thanks.

Improve this question
$\endgroup$
3
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – Michael E2 Feb 28 at 20:05
  • 1
    $\begingroup$ Why not to put $t_0=0$? $\endgroup$ – Alex Trounev Feb 28 at 21:16
  • $\begingroup$ Yes t_0 can be 0. $\endgroup$ – Moe Mar 1 at 6:06
4
$\begingroup$

I think this should do it. First define the piecewise function:

g[u_List, tmax_] := Block[{t},
  Function[t,
   Evaluate @ Dot[u,
     Boole[#1 <= t < #2] & @@@
      Partition[Subdivide[0, tmax, Length[u]], 2, 1]
     ]
   ]
  ]

g is a function that returns a function. For example:

g[{1, 0, 1}, 1]
%[0.8]

Next, define the equations. I added an extra equation that keeps track of the integral using the variable int:

Clear[x, y, u, t, int, uvec, npts]
npts = 5;
tmax = 1;
uvec = Array[u, npts];
eqs = {
  x'[t] == 2 x[t] - x[t] y[t] - 3 x[t] g[uvec, tmax][t],
  y'[t] == y[t] - 2 x[t] y[t] - x[t] g[uvec, tmax][t],
  int'[t] == (x[t] - 1)^2 + (y[t] - 1)^2,
  x[0] == 1, y[0] == 1/2, int[0] == 0
  }

Solve using ParametricNDSolveValue:

sol = ParametricNDSolveValue[
  eqs,
  {int[tmax], Function[t, {x[t], y[t]}]},
  {t, 0, tmax},
  uvec
  ]

Test if it works for a typical u vector:

sol[1, 1, 0, 0, 1]

{0.7955, Function[t$, ...]}

Minimize the solution over the given constraints:

min = NMinimize[
  {
   Indexed[sol @@ uvec, 1],
   uvec \[Element] Integers && And @@ Map[# == 0 || # == 1 &, uvec]
   },
  uvec
 ]

{0.7955, {u[1] -> 1, u[2] -> 1, u[3] -> 0, u[4] -> 0, u[5] -> 1}}

Plot the solution:

ParametricPlot[
 Evaluate[Last[sol @@ uvec /. Last[min]][t]],
 {t, 0, tmax}
]
Plot[
  Evaluate[Last[sol @@ uvec /. Last[min]][t]], {t, 0, tmax}, 
  PlotLabel -> {uvec /. min[[2]]}
]

enter image description here

enter image description here

Solution for npts = 8:

enter image description here

Edit

It should be noted that this approach scales very badly because the discrete nature of the u[i] doesn't allow NMinimize to use any sort of gradient-based method. This means that it just has to brute-force all 2^npts possibilities.

Instead, consider something like:

min = FindMinimum[
  Evaluate @ {
     Indexed[sol @@ uvec, 1],
     uvec \[Element] Reals && And @@ Map[0 <= # <= 1 &, uvec]
  },
  Evaluate @ Thread[{uvec, 1/2}]
]
Improve this answer
$\endgroup$
10
  • 1
    $\begingroup$ (+1) this is a good answer. Also for npts=8 it is 1:1 to solution I got with different algorithm without NDSolve and NMinimize. Can you add picture Plot[Evaluate[Last[sol @@ uvec /. Last[min]][t]], {t, 0, tmax}, PlotLabel -> {uvec /. min[[2]]}] ? $\endgroup$ – Alex Trounev Mar 1 at 13:52
  • $\begingroup$ @AlexTrounev Pictures has been added. It looks like the function g[t] can stabilize the x coordinate reasonably well, but the y coordinate just runs away. $\endgroup$ – Sjoerd Smit Mar 1 at 14:08
  • $\begingroup$ Thank you! I mean picture with npts=8 to compare then with my code. Now I try to get solution with npts=16 to test you code. $\endgroup$ – Alex Trounev Mar 1 at 14:33
  • $\begingroup$ @AlexTrounev Added pictures for 8 control points. $\endgroup$ – Sjoerd Smit Mar 1 at 16:04
  • $\begingroup$ Thank you! I have added code to compare. Can you run your code with npts=16? $\endgroup$ – Alex Trounev Mar 1 at 18:26
3
$\begingroup$

This problem can be solved with collocation method using Bernoulli wavelets. We don't need to run NDSolve and NMinimize so it is alternative to that proposed by Sjoerd Smit. Nevertheless we can compare final result for npts=8. First we define wavelets, functions and derivatives as follows

n = 2;


M = Sum[1, {j, 0, n, 1}, {i, 0, 2^j - 1, 1}] + 1; U = Array[u, {M}];
dx = 1/M; A = 0; xl = Table[A + l*dx, {l, 0, M}]; tcol = 
 xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, M + 1}]; 
psi1[x_] := Piecewise[{{BernoulliB[2, x], 0 <= x < 1}, {0, True}}];
psi2[x_] := Piecewise[{{BernoulliB[1, x], 0 <= x < 1}, {0, True}}];
psi1jk[x_, j_, k_] := psi1[j*x - k];
psi2jk[x_, j_, k_] := psi2[j*x - k];
psijk = Compile[{{x, _Real}, {j, _Integer}, {k, _Integer}}, (psi1jk[x,
        j, k] + psi2jk[x, j, k])/2];
(*intjk=Integrate[psijk[x,j,k],x,Assumptions\[Rule]{j>0,k>0}]*)
psijk1 = Compile[{{x, _Real}, {j, _Integer}, {k, _Integer}}, 
   Piecewise[{{(-k + k^2)/(2*j), 
      j > 0 && k == 0 && 
       1/j - x < 0}, {(1/6)*(-x + 3*k^2*x - 3*j*k*x^2 + j^2*x^3), 
      j > 0 && k == 0 && x > 0 && 
       1/j - x >= 
        0}, {(k - k^3 - j*x + 3*j*k^2*x - 3*j^2*k*x^2 + j^3*x^3)/(6*
         j), j > 0 && k > 0 && k/j - x < 0 && 1/j + k/j - x >= 0}}, 
    0]];
Psi[x_] := 
  Join[{1}, 
   Flatten[Table[psijk[x, 2^j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]];
int1[x_] := 
  Join[{x}, 
   Flatten[Table[
     psijk1[x, 2^j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]];
var1 = Join[{a0}, 
  Flatten[Table[a[j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]]; var2 = 
 Join[{b0}, 
  Flatten[Table[b[j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]]; var3 = 
 Join[{c0}, Flatten[Table[c[j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]];
z1[t_] := var3.Psi[t]; z[t_] := var3.int1[t] + c1;
y1[t_] := var1.Psi[t]; y[t_] := var1.int1[t] + a1; 
x1[t_] := var2.Psi[t]; x[t_] := var2.int1[t] + b1;

Second, we define model to be optimized

varM = Join[{a1, b1, c1}, var1, var2, var3]; tf = 1; 
j = 0; 
AbsoluteTiming[Do[j = j + 1; eq = Flatten[Table[{2*x[xcol[[i]]] - 3*u[i]*x[xcol[[i]]] - x1[xcol[[i]]]/tf - x[xcol[[i]]]*y[xcol[[i]]] == 0, 
        (-u[i])*x[xcol[[i]]] + y[xcol[[i]]] - 2*x[xcol[[i]]]*y[xcol[[i]]] - y1[xcol[[i]]]/tf == 0, 
        -z1[xcol[[i]]]/tf + (x[xcol[[i]]] - 1)^2 + (y[xcol[[i]]] - 1)^2 == 0}, {i, Length[xcol]}]]; 
    cons = Join[eq, {x[0] == 1, y[0] == 1/2, z[0] == 0}]; sol[j] = Quiet[FindRoot[cons, Table[{varM[[i]], 1/10}, {i, Length[varM]}]]]; 
    s[j] = z[1] /. sol[j]; uu[j] = U; , {u[1], 0, 1, 1}, {u[2], 0, 1, 1}, {u[3], 0, 1, 1}, {u[4], 0, 1, 1}, {u[5], 0, 1, 1}, {u[6], 0, 1, 1}, 
   {u[7], 0, 1, 1}, {u[8], 0, 1, 1}]]

Third, we can fined minimum and plot solution as follows

{km, sm} = MinimalBy[Table[{k, s[k]}, {k, 1, j}], Last][[1]]
(*{227, 0.786878}*)
lst1 = Join[{{0, 1}}, 
  Table[{xcol[[i]], x[xcol[[i]]] /. sol[km]}, {i, M}]]; lst2 = 
 Join[{{0, 1/2}}, Table[{xcol[[i]], y[xcol[[i]]] /. sol[km]}, {i, M}]];

ListLinePlot[{lst1, lst2}, PlotLabel -> uu[km], Frame -> True, 
 Axes -> False, PlotLegends -> {"x", "y"}]

Figure 1

So far so good, we have practically same minimum value 0.786878 compare to solution by Sjoerd Smit of 0.787574. Now we put n=3 in my code that corresponds to npts = 16. In this case model to be optimized is

    varM = Join[{a1, b1, c1}, var1, var2, var3]; tf = 1; 
j = 0; 
AbsoluteTiming[Do[j = j + 1; eq = Flatten[Table[{2*x[xcol[[i]]] - 3*u[i]*x[xcol[[i]]] - x1[xcol[[i]]]/tf - x[xcol[[i]]]*y[xcol[[i]]] == 0, 
        (-u[i])*x[xcol[[i]]] + y[xcol[[i]]] - 2*x[xcol[[i]]]*y[xcol[[i]]] - y1[xcol[[i]]]/tf == 0, 
        -z1[xcol[[i]]]/tf + (x[xcol[[i]]] - 1)^2 + (y[xcol[[i]]] - 1)^2 == 0}, {i, Length[xcol]}]]; 
    cons = Join[eq, {x[0] == 1, y[0] == 1/2, z[0] == 0}]; sol[j] = Quiet[FindRoot[cons, Table[{varM[[i]], 1/10}, {i, Length[varM]}]]]; 
    s[j] = z[1] /. sol[j]; uu[j] = U; , {u[1], 0, 1, 1}, {u[2], 0, 1, 1}, {u[3], 0, 1, 1}, {u[4], 0, 1, 1}, {u[5], 0, 1, 1}, {u[6], 0, 1, 1}, 
   {u[7], 0, 1, 1}, {u[8], 0, 1, 1}, {u[9], 0, 1, 1}, {u[10], 0, 1, 1}, {u[11], 0, 1, 1}, {u[12], 0, 1, 1}, {u[13], 0, 1, 1}, 
   {u[14], 0, 1, 1}, {u[15], 0, 1, 1}, {u[16], 0, 1, 1}]]

With my code we got final result in 26 min at j=65536. New result is

 {km, sm} = MinimalBy[Table[{k, s[k]}, {k, 1, j}], Last][[1]]

(*Out[]= {63763, 0.782377} *)

This result not so differ from above while picture not looks similar

Figure 2

This code can be optimize with Compile and ParallelDo as it shown here. We also can compare last picture with NMimimize (it finished after 8 h 27 min running on my laptop). Final result looks like

min = NMinimize[{Indexed[sol @@ uvec, 1], 
        uvec \[Element] Integers && And @@ Map[# == 0 || # == 1 &, uvec]},
        uvec] // AbsoluteTiming
    
Out[]= {30454.7, {0.783045, {u[1] -> 1, u[2] -> 1, u[3] -> 1, 
   u[4] -> 1, u[5] -> 1, u[6] -> 0, u[7] -> 0, u[8] -> 0, u[9] -> 1, 
   u[10] -> 0, u[11] -> 0, u[12] -> 1, u[13] -> 0, u[14] -> 0, 
   u[15] -> 1, u[16] -> 0}}}

Now we have difference in vector U as well, while minimum value with my code and with NMinimize - 0.782377 and 0.783045 consequently. Pictures also are differ but it can be explained by different interpolation method with wavelets and NDSolve. Figure 4

Improve this answer
$\endgroup$
9
  • $\begingroup$ Thanks to all for valuable insights. My ultimate goal is to see if Mathematica NMinimize or FindMinimum can deal with a large scale problem where I ultimately want to consider the discritization step from t0 to tf to be fine enough such as dt=10^(-3) and this will produce 2^(10^3) possibilities. That is HUGE !!! $\endgroup$ – Moe Mar 3 at 7:36
  • $\begingroup$ Also, I want to comment on Bernoulli wavelets approach. Why specifically such basis functions? Do they have specific approximation property over other basis functions? $\endgroup$ – Moe Mar 3 at 7:37
  • $\begingroup$ @Moe With wavelets approach it is possible to make small step size, but with NMinimize it is not possible to solve this problem directly since even for dt=1/16 it takes 8 h 27 min while with wavelets it takes 26 min without optimization. Also it should be algorithm different from that shown above. I am testing Bernoulli wavelets for different problems like this one mathematica.stackexchange.com/questions/226421/… $\endgroup$ – Alex Trounev Mar 3 at 10:37
  • $\begingroup$ Great! Can these wavelets (or the approach of using them) handle a fine step size say h=0.01 with say t0=0 and tmax=12.? As I mentioned in my previous comment, my ultimate goal to consider the problem of having number of control points as a parameter as well as the final time tmax, in addition to u[i]. So I would have to minimize an cost function with respect to control points t1, t2, t3, ...,tN, (where N is also a parameter) and the u[i], i=1,...,N $\endgroup$ – Moe Mar 4 at 6:52
  • $\begingroup$ @Moe Problem with unknown $N$ is different from that we discussed here. Also it can be solved with wavelets and with algorithm we need to develop. This new algorithm should by supported with some theorem about convergence of such solution. $\endgroup$ – Alex Trounev Mar 4 at 12:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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