# Local optimization of parameters in a Coupled Non-Linear System

I have a system of differential equations :

$\frac{dX}{dt} = (0.8+0.4c_{1})X-(0.007+0.006c_{2})XY$ and

$\frac{dY}{dt} = (0.048+0.004c_{3})XY-(0.4+0.2c_{4})Y$

where $0\le c_{i}\le1$.

Varying $c_{i}$ gives me an envelope of solutions for both $X(t)$ and $Y(t)$. My goal is to find the combination of $c_{i}$ which gives the upper and lower limits of my solution.

Taking $c_{i}$ as discrete points between 0 and 1 with step size 0.1, I used the following code to obtain a solution, and managed to plot it.

sol = ParametricNDSolve[{x'[t] == (0.8 + 0.4 c1) x[t] - (0.007 + 0.006 c2) x[t] y[t],y'[t] == (0.0048 + 0.0004 c3) x[t] y[t] - (0.4 + 0.2 c4) y[t], x == 200, y == 80}, {x, y}, {t, 0, 50}, {c1, c2, c3, c4}]
xt = Flatten[Evaluate[Table[x[c1, c2, c3, c4][t] /. sol, {c1, 0, 1, 0.1}, {c2, 0, 1,0.1}, {c3, 0, 1, 0.1}, {c4, 0, 1, 0.1}]]];
mat = Table[Table[Evaluate[xt[[a]]], {t, 0, 25, 0.01}], {a,Range[Length[xt]]}];
l = Min /@ Transpose[mat];
m = Max /@ Transpose[mat];
ind = Range[0, 25, 0.01];
lupd = Transpose[{ind, l}];
mupd = Transpose[{ind, m}];
ListLinePlot[{lupd, mupd},GridLines -> Automatic, ImageSize -> {600, 300}]


The plot for $X(t)$ came as The plot for $Y(t)$ came as Edit 3

It is not possible to calculate the equation of the upper curve of the envelope obtained by varying $c_{i}$s. Hence, I reworked my question.

Consider the following code fragment

Plot[Evaluate@({x[1, 0, 0, 1][t], x[0, 1, 1, 1][t], x[0, 0, 0, 0][t],
x[0.8, 0, 0, 0.2][t]} /. sol), {t, 0, 25}]


This gives 4 different trajectories for $X(t)$. Consider the line $t=12$. As it cuts the enevelope, we get a minimum point and a maximum point.

I would like the combination of $c_{i}$ that give this maximum and minimum points for $t=1,2,3,...,25$.

In my solution I used Max and Min to extract these maximum and minimum points and then plotted them. I could not get the combination of $c_{i}$ for each of these points.

• Any optimization does require that you construct a merit function; that is, what exactly are you optimizing? What should be "big", and what should be "small"? – J. M.'s discontentment Jan 24 '17 at 16:53
• The solution of the coupled equation. $X(t)$ and $Y(t)$ in this case. – Aritra Saha Jan 24 '17 at 16:59
• How can a time-varying oscillating solution be "small"? Isn't this the lower bounds (blue) curves that you are showing on the plots? – gpap Jan 24 '17 at 17:04
• In short, you effectively want to maximize $\|X(t)\|_\infty$ and $\|Y(t)\|_\infty$ in the "upper envelope" case? – J. M.'s discontentment Jan 24 '17 at 17:14
• It is possible that the $c_i$ that maximize $\|X(t)\|_\infty$ can be different from the $c_i$ that maximize $\|Y(t)\|_\infty$; what do you propose to do in that case? – J. M.'s discontentment Jan 24 '17 at 17:20

As I pointed out in the comments, the $c_i$ do not retain the frequency of the solutions constant so the solutions cross each other and there isn't one collection of $c_i$ that satisfies the "maximum envelope" criterion. For different integer times it's an easy application of NMaximize:

Table[NMaximize[{x[c1, c2, c3, c4][t] /. sol, 0 <= c1 <= 1, 0 <= c2 <= 1, 0 <= c3 <= 1, 0 <= c4 <= 1}, {c1, c2, c3, c4}], {t, 1, 25}]


which produces a solution (but doesn't always converge properly so maybe some playing around with the parameters of NMaximize is necessary here):

{{311.231, {c1 -> 1., c2 -> 0., c3 -> 0.,
c4 -> 1.}}, {195.645, {c1 -> 0.851284, c2 -> 0., c3 -> 0.,
c4 -> 1.}}, {101.746, {c1 -> 0., c2 -> 0., c3 -> 0.,
c4 -> 1.}}, {72.4078, {c1 -> 0.554206, c2 -> 1., c3 -> 0.,
c4 -> 1.}}, {86.7196, {c1 -> 0.516399, c2 -> 1., c3 -> 0.,
c4 -> 1.}}, {119.434, {c1 -> 0.783024, c2 -> 1., c3 -> 0.,
c4 -> 1.}}, {174.87, {c1 -> 1., c2 -> 1., c3 -> 0.,
c4 -> 1.}}, {225.381, {c1 -> 1., c2 -> 0.393763, c3 -> 0.,
c4 -> 1.}}, {312.48, {c1 -> 1., c2 -> 0., c3 -> 0.,
c4 -> 1.}}, {298.507, {c1 -> 1., c2 -> 0., c3 -> 0.,
c4 -> 0.483703}}, {289.22, {c1 -> 1., c2 -> 0., c3 -> 0.,
c4 -> 0.0869823}}, {242.738, {c1 -> 0.536497, c2 -> 0., c3 -> 0.,
c4 -> 0.}}, {206.276, {c1 -> 0., c2 -> 4.11611*10^-9,
c3 -> 8.75751*10^-6, c4 -> 0.}}, {143.553, {c1 -> 1., c2 -> 1.,
c3 -> 0., c4 -> 1.}}, {196.936, {c1 -> 1., c2 -> 1., c3 -> 0.,
c4 -> 1.}}, {245.969, {c1 -> 1., c2 -> 0.331846, c3 -> 0.,
c4 -> 1.}}, {313.172, {c1 -> 1., c2 -> 0., c3 -> 0.,
c4 -> 1.}}, {304.824, {c1 -> 1., c2 -> 0., c3 -> 0.,
c4 -> 0.715149}}, {298.267, {c1 -> 1., c2 -> 0., c3 -> 0.,
c4 -> 0.4678}}, {293.266, {c1 -> 1., c2 -> 0.0000225295,
c3 -> 2.52908*10^-7, c4 -> 0.262797}}, {281.43, {c1 -> 0.99607,
c2 -> 0.00393228, c3 -> 0.739327,
c4 -> 0.108152}}, {270.711, {c1 -> 0.939643, c2 -> 0.000403377,
c3 -> 0.999011, c4 -> 0.000087475}}, {211.763, {c1 -> 1.,
c2 -> 0.757642, c3 -> 0., c4 -> 1.}}, {261.286, {c1 -> 1.,
c2 -> 0.265911, c3 -> 0., c4 -> 1.}}, {313.282, {c1 -> 1.,
c2 -> 0., c3 -> 0., c4 -> 1.}}}

• What is the runtime for this code segment? Also, how do I tabulate this result? (I have to include of a table of $c_{i}$s for different $t$ in my report) – Aritra Saha Jan 27 '17 at 15:05
• This took about 10sec on my computer. For making a table I think this is relevant – gpap Jan 27 '17 at 15:08
• I will go through the link. How do you think the increase in the number of parameters will affect the runtime? – Aritra Saha Jan 27 '17 at 15:12
• It will increase it :) By how much depends on how many more parameters, how long a time etc. – gpap Jan 27 '17 at 15:15
• I think that is a separate question in itself but you are essentially solving a nonlinear ode system $10^4$ times to obtain the solutions for each parameter – gpap Jan 27 '17 at 15:32