4
$\begingroup$

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"]
Improve this question
$\endgroup$
1

1 Answer 1

Reset to default
2
$\begingroup$

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.

Improve this answer
$\endgroup$

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.