Tell me more ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.
up vote 2 down vote favorite

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"]
share|improve this question
1  
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

1 Answer

up vote 1 down vote

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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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