# Minimization using NMinimize

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.

• 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! Commented Feb 28, 2021 at 20:05
• Why not to put $t_0=0$? Commented Feb 28, 2021 at 21:16
• Yes t_0 can be 0.
– Moe
Commented Mar 1, 2021 at 6:06

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]