# Problems with FindMinimum for Constrained Optimization

I am trying to define a simple optimization problem, arising from the discretization of a simple Optimal Control problem.

The problem that results is essentially linear, with bilinearities only appearing in some equality constraints. The problem is not big and should be solvable by FindMinimum.

However I get the following message:

Power::infy: Infinite expression 1/0. encountered. >>
Min::nord: Invalid comparison with ComplexInfinity attempted. >>
FindMinimum::conv: Interior point method fails to converge. >>


I note that the constraints do not contain divisions at all, and the objective is totally linear.

The complete model is given below:

Simple car problem
\begin{align*}\frac{dv}{dt} &= u\\ \frac{dx}{dt} &= v\end{align*}
$x(0)=0,\quad v(0)=0,\quad x(t_f)=300,\quad v(t_f)=0, \qquad \text{minimize } t_f$.
$-2 \leq u\leq +1$.

In[51]:= (* try to define the car problem directly below *)
(* v is velocity and x is distance, h is the stepsizes, u are the controls *)

nsteps = 10;

vlist = Table[v[i], {i, 0, nsteps}]

Out[52]= {v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8], v[9], v[10]}

In[53]:= xlist = Table[x[i], {i, 0, nsteps}]

Out[53]= {x[0], x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10]}

In[54]:= hlist = Table[h[i], {i, 1, nsteps}]

Out[54]= {h[1], h[2], h[3], h[4], h[5], h[6], h[7], h[8], h[9], h[10]}

In[55]:= ulist = Table[u[i], {i, 1, nsteps}]

Out[55]= {u[1], u[2], u[3], u[4], u[5], u[6], u[7], u[8], u[9], u[10]}

In[56]:= varslist = Join[vlist, xlist, hlist, ulist]

Out[56]= {v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8], v[9], v[10], x[0],
x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10], h[1], h[2],
h[3], h[4], h[5], h[6], h[7], h[8], h[9], h[10], u[1], u[2], u[3], u[4],
u[5], u[6], u[7], u[8], u[9], u[10]}

In[57]:= obj0 = Sum[h[i], {i, 1, nsteps}]

Out[57]= h[1] + h[2] + h[3] + h[4] + h[5] + h[6] + h[7] + h[8] + h[9] + h[10]

In[58]:= odes = Flatten[
Table[{x[i] - x[i - 1] == v[i]*h[i], v[i] - v[i - 1] == h[i]*u[i]}, {i, 1,
nsteps}]]

Out[58]= {-x[0] + x[1] == h[1] v[1], -v[0] + v[1] == h[1] u[1], -x[1] + x[2] ==
h[2] v[2], -v[1] + v[2] == h[2] u[2], -x[2] + x[3] ==
h[3] v[3], -v[2] + v[3] == h[3] u[3], -x[3] + x[4] ==
h[4] v[4], -v[3] + v[4] == h[4] u[4], -x[4] + x[5] ==
h[5] v[5], -v[4] + v[5] == h[5] u[5], -x[5] + x[6] ==
h[6] v[6], -v[5] + v[6] == h[6] u[6], -x[6] + x[7] ==
h[7] v[7], -v[6] + v[7] == h[7] u[7], -x[7] + x[8] ==
h[8] v[8], -v[7] + v[8] == h[8] u[8], -x[8] + x[9] ==
h[9] v[9], -v[8] + v[9] == h[9] u[9], -x[9] + x[10] ==
h[10] v[10], -v[9] + v[10] == h[10] u[10]}

In[59]:= ubounds = Flatten[Table[{u[i] <= 1.0, u[i] >= -2.0}, {i, 1, nsteps}]]

Out[59]= {u[1] <= 1., u[1] >= -2., u[2] <= 1., u[2] >= -2., u[3] <= 1., u[3] >= -2.,
u[4] <= 1., u[4] >= -2., u[5] <= 1., u[5] >= -2., u[6] <= 1., u[6] >= -2.,
u[7] <= 1., u[7] >= -2., u[8] <= 1., u[8] >= -2., u[9] <= 1., u[9] >= -2.,
u[10] <= 1., u[10] >= -2.}

In[60]:= hbounds = Flatten[Table[{h[i] <= 100.0, h[i] >= 0.01}, {i, 1, nsteps}]]

Out[60]= {h[1] <= 100., h[1] >= 0.01, h[2] <= 100., h[2] >= 0.01, h[3] <= 100.,
h[3] >= 0.01, h[4] <= 100., h[4] >= 0.01, h[5] <= 100., h[5] >= 0.01,
h[6] <= 100., h[6] >= 0.01, h[7] <= 100., h[7] >= 0.01, h[8] <= 100.,
h[8] >= 0.01, h[9] <= 100., h[9] >= 0.01, h[10] <= 100., h[10] >= 0.01}

In[61]:= boundary = {x[0] == 0.0, v[0] == 0.0, v[nsteps] == 0.0, x[nsteps] == 300.0}

Out[61]= {x[0] == 0., v[0] == 0., v[10] == 0., x[10] == 300.}

In[62]:= constr0 = Join[odes, ubounds, hbounds, boundary]

Out[62]= {-x[0] + x[1] == h[1] v[1], -v[0] + v[1] == h[1] u[1], -x[1] + x[2] ==
h[2] v[2], -v[1] + v[2] == h[2] u[2], -x[2] + x[3] ==
h[3] v[3], -v[2] + v[3] == h[3] u[3], -x[3] + x[4] ==
h[4] v[4], -v[3] + v[4] == h[4] u[4], -x[4] + x[5] ==
h[5] v[5], -v[4] + v[5] == h[5] u[5], -x[5] + x[6] ==
h[6] v[6], -v[5] + v[6] == h[6] u[6], -x[6] + x[7] ==
h[7] v[7], -v[6] + v[7] == h[7] u[7], -x[7] + x[8] ==
h[8] v[8], -v[7] + v[8] == h[8] u[8], -x[8] + x[9] ==
h[9] v[9], -v[8] + v[9] == h[9] u[9], -x[9] + x[10] ==
h[10] v[10], -v[9] + v[10] == h[10] u[10], u[1] <= 1., u[1] >= -2.,
u[2] <= 1., u[2] >= -2., u[3] <= 1., u[3] >= -2., u[4] <= 1., u[4] >= -2.,
u[5] <= 1., u[5] >= -2., u[6] <= 1., u[6] >= -2., u[7] <= 1., u[7] >= -2.,
u[8] <= 1., u[8] >= -2., u[9] <= 1., u[9] >= -2., u[10] <= 1., u[10] >= -2.,
h[1] <= 100., h[1] >= 0.01, h[2] <= 100., h[2] >= 0.01, h[3] <= 100.,
h[3] >= 0.01, h[4] <= 100., h[4] >= 0.01, h[5] <= 100., h[5] >= 0.01,
h[6] <= 100., h[6] >= 0.01, h[7] <= 100., h[7] >= 0.01, h[8] <= 100.,
h[8] >= 0.01, h[9] <= 100., h[9] >= 0.01, h[10] <= 100., h[10] >= 0.01,
x[0] == 0., v[0] == 0., v[10] == 0., x[10] == 300.}

In[63]:= nvars0 = Dimensions[varslist][[1]]

Out[63]= 42

In[64]:= vars0 = Transpose[{varslist, Table[Random[Real, {0., 1.}], {nvars0}]}]

Out[64]= {{v[0], 0.994273}, {v[1], 0.132455}, {v[2], 0.7492}, {v[3], 0.422551}, {v[4],
0.3127}, {v[5], 0.734046}, {v[6], 0.386862}, {v[7], 0.95813}, {v[8],
0.319134}, {v[9], 0.0195711}, {v[10], 0.547824}, {x[0], 0.126578}, {x[1],
0.91368}, {x[2], 0.385664}, {x[3], 0.195437}, {x[4], 0.589594}, {x[5],
0.369296}, {x[6], 0.437719}, {x[7], 0.205286}, {x[8], 0.210231}, {x[9],
0.145847}, {x[10], 0.225091}, {h[1], 0.572089}, {h[2], 0.681713}, {h[3],
0.151574}, {h[4], 0.0926365}, {h[5], 0.822889}, {h[6], 0.259161}, {h[7],
0.838874}, {h[8], 0.35859}, {h[9], 0.436027}, {h[10], 0.301031}, {u[1],
0.51974}, {u[2], 0.339019}, {u[3], 0.888204}, {u[4], 0.174453}, {u[5],
0.60606}, {u[6], 0.953355}, {u[7], 0.692767}, {u[8], 0.584859}, {u[9],
0.236764}, {u[10], 0.515636}}

In[65]:= constr0 = Apply[And, constr0]

Out[65]= -x[0] + x[1] == h[1] v[1] && -v[0] + v[1] == h[1] u[1] && -x[1] + x[2] ==
h[2] v[2] && -v[1] + v[2] == h[2] u[2] && -x[2] + x[3] ==
h[3] v[3] && -v[2] + v[3] == h[3] u[3] && -x[3] + x[4] ==
h[4] v[4] && -v[3] + v[4] == h[4] u[4] && -x[4] + x[5] ==
h[5] v[5] && -v[4] + v[5] == h[5] u[5] && -x[5] + x[6] ==
h[6] v[6] && -v[5] + v[6] == h[6] u[6] && -x[6] + x[7] ==
h[7] v[7] && -v[6] + v[7] == h[7] u[7] && -x[7] + x[8] ==
h[8] v[8] && -v[7] + v[8] == h[8] u[8] && -x[8] + x[9] ==
h[9] v[9] && -v[8] + v[9] == h[9] u[9] && -x[9] + x[10] ==
h[10] v[10] && -v[9] + v[10] == h[10] u[10] && u[1] <= 1. && u[1] >= -2. &&
u[2] <= 1. && u[2] >= -2. && u[3] <= 1. && u[3] >= -2. && u[4] <= 1. &&
u[4] >= -2. && u[5] <= 1. && u[5] >= -2. && u[6] <= 1. && u[6] >= -2. &&
u[7] <= 1. && u[7] >= -2. && u[8] <= 1. && u[8] >= -2. && u[9] <= 1. &&
u[9] >= -2. && u[10] <= 1. && u[10] >= -2. && h[1] <= 100. && h[1] >= 0.01 &&
h[2] <= 100. && h[2] >= 0.01 && h[3] <= 100. && h[3] >= 0.01 &&
h[4] <= 100. && h[4] >= 0.01 && h[5] <= 100. && h[5] >= 0.01 &&
h[6] <= 100. && h[6] >= 0.01 && h[7] <= 100. && h[7] >= 0.01 &&
h[8] <= 100. && h[8] >= 0.01 && h[9] <= 100. && h[9] >= 0.01 &&
h[10] <= 100. && h[10] >= 0.01 && x[0] == 0. && v[0] == 0. && v[10] == 0. &&
x[10] == 300.

In[66]:= FindMinimum[
{obj0, constr0}, vars0, Method -> "InteriorPoint"
]

During evaluation of In[66]:= Power::infy: Infinite expression 1/0. encountered. >>

During evaluation of In[66]:= Min::nord: Invalid comparison with ComplexInfinity attempted. >>

During evaluation of In[66]:= FindMinimum::conv: Interior point method fails to converge. >>

Out[66]= FindMinimum[{obj0, constr0}, vars0, Method -> "InteriorPoint"]

-
Not sure but it seems that your problem is related to a similar one which was answered here : mathematica.stackexchange.com/questions/10362/…. – Artes Sep 13 '12 at 0:38

I'm not sure, but this might help a little:

Your system can be written as

$$\frac{d^2 x}{d t^2} = u(x),$$

and with $x(0) = 0$, $x(t_f) = 300$ it has as solution [1]

$$x(t) = \frac{300 t}{t_f} + \frac{t-t_f}{t_f}\int_0^t s u(s) ds + \frac{t}{t_f}\int_t^{t_f} (s-t_f)u(s) ds,$$

as one can verify with MMA

x[t_, tf_] := 300 t/tf + (t - tf)/tf Integrate[s u[s], {s, 0, t}] +
t/tf Integrate[(s - tf) u[s], {s, t, tf}]

Simplify@{x[0, tf], x[tf, tf], D[x[t, tf], {t, 2}]}

(* {0, 300, u[t]} *)


Then

$$v(t_f) = \frac{1}{t_f}\left(300 + \int_0^{t_f} s u(s) ds \right)$$

And hence, the problem can be rephrased as

Find the minimum of $f(t_f) = t_f$ subject to the constraint $$0 = \left(300 + \int_0^{t_f} s u(s) ds \right)$$

Now, I don't fully understand what you are doing and I don't know much on control theory, but maybe you can try to solve the problem using a trapezoidal approximation along with your random variables.

[1] Friedman B, Principles and Techniques of Applied Mathematics, Wiley & Sons, 1956.

-