<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<style type="text/css" style="display:none;"><!-- P {margin-top:0;margin-bottom:0;} --></style>
</head>
<body dir="ltr">
<div id="divtagdefaultwrapper" style="font-size:12pt;color:#000000;font-family:Calibri,Arial,Helvetica,sans-serif;">
<p>Dear all,</p>
<p><br>
</p>
<p>I am having a problem using Ipopt: After having modelling a simple dynamic system based on newton rules and discretize it, I am getting a maximun number of iterations exceeded problem. I am not sure if the problem is due to the discretization method (a leapfrog
 integration method), because using a simple Euler method (look for "##" in the code below) the system convergency is got in a few iterations.</p>
<p><br>
</p>
<p>Does someone understand where the problem is? It is worth the extra precission of having into acount the aceleration for calculating the space or the euler method would be good enought for prototyping?</p>
<p><br>
</p>
<p>Please, I copy the AMPL code I am using. Take into account that the low number of steps is due to the restricted version of AMPL. Anyway, using a bigger number in Neos server does not solve the problem.</p>
<p><br>
</p>
<p>I am not sure if this mailling list is for this kind of questions. If it is not I would be sorry, please, could someone tell me about any forum of Ipopt in particular and optimal control in general?</p>
<p><br>
</p>
<p>Thanks in advance,</p>
<p>Pablo</p>
<p><br>
</p>
<p><br>
</p>
<p></p>
<div># We have a train hat is at position x=0 at time t=0, and we can</div>
<div># control the acceleration, a, (which includes braking as well).</div>
<div>#</div>
<div># N:  Number of iterations.</div>
<div># L:  Path length.</div>
<div># T:  Travel time initial guess.</div>
<div># aU: Upper bound on acceleration.</div>
<div># aL: Lower bound on acceleration.</div>
<div># jU: Upper bound on jerk.</div>
<div># jL: Lower bound on jerk.</div>
<div># Ra: Drag coeficient (independent).</div>
<div># Rb: Drag coeficient (v dependant).</div>
<div># Rc: Drag coeficient (v^2 dependant).</div>
<div>#</div>
<div># The overall optimal control problem is therefore:</div>
<div>#</div>
<div># min  tf</div>
<div># s.t. dx/dt(t) = v(t)                               for t in [0,tf]</div>
<div>#      dv/dt(t) = a(t) - Ra - Rb*v(t) - Rc*v(t)^2    for t in [0,tf]</div>
<div>#      aL   <= dv/dt(t)   <=  aU                     for t in [0,tf]</div>
<div>#      jL   <= d2v/dt2(t) <=  jU                     for t in [0,tf]</div>
<div>#      v(t) <= Vmax(x)                               for t in [0,tf]</div>
<div>#</div>
<div>#      x(0)  = 0; v(0)  = 0</div>
<div>#      x(tf) = L; v(tf) = 0</div>
<div>#</div>
<div># Leapfrog discretization.</div>
<div># With this, we can state the overall optimization problem below:</div>
<div><br>
</div>
<div>option solver ipopt;</div>
<div><br>
</div>
<div># Number of discretization intervals</div>
<div>param N := 1000;</div>
<div><br>
</div>
<div># Final position</div>
<div>param L := 1000;</div>
<div><br>
</div>
<div># Final time</div>
<div>param T := 50;</div>
<div><br>
</div>
<div># Upper and lower bound on acceleration</div>
<div>param aU := 1.15;</div>
<div>param aL := -1.2;</div>
<div><br>
</div>
<div># Upper and lower bound on jerk</div>
<div>param jU := 0.8;</div>
<div>param jL := -0.8;</div>
<div><br>
</div>
<div># Drag factors</div>
<div>param Ra := 5.71428/322;</div>
<div>param Rb := 0.11688/322;</div>
<div>param Rc := 0.00842/322;</div>
<div><br>
</div>
<div># Size of discretization intervals</div>
<div>var tf >= 0 := T;</div>
<div><br>
</div>
<div>var h = tf/N;</div>
<div><br>
</div>
<div># Velocity</div>
<div>var x{i in 0..N} >= 0, := i*(L/N);</div>
<div><br>
</div>
<div># Time</div>
<div>var v{i in 0..N} >= 0, := L/T;</div>
<div><br>
</div>
<div># Acceleration</div>
<div>var a{i in 0..N} := 0;</div>
<div><br>
</div>
<div># Objective function</div>
<div>minimize time:</div>
<div>    tf;</div>
<div><br>
</div>
<div>subject to Amax {i in 0..N}:</div>
<div>    a[i] <= if (v[i] <= 11) then (343/322) else (343/322)*(11/v[i]);</div>
<div>subject to Dmax {i in 0..N}:</div>
<div>    a[i] >= if (v[i] <= 19) then -(327/322) else -(327/322)*(19/v[i]);</div>
<div>    </div>
<div># Differential equation for velocity</div>
<div>subject to dx {i in 0..N-1}:</div>
<div>    (x[i+1]-x[i])/h = v[i] + 0.5*(a[i] - Ra - Rb*v[i] - Rc*v[i]^2)*h;</div>
<div>    ##(x[i+1]-x[i])/h = v[i];</div>
<div>    </div>
<div># Differential equation for time</div>
<div>subject to dv {i in 0..N-1}:</div>
<div>    (v[i+1]-v[i])/h = a[i] - Ra - Rb*v[i] - Rc*v[i]^2;</div>
<div><br>
</div>
<div># Confort constraints</div>
<div>subject to acel {i in 0..N-1}:</div>
<div>    aL <= (v[i+1]-v[i])/h <= aU;</div>
<div>    </div>
<div>subject to jerk {i in 1..N-1}:</div>
<div>    jL <= (v[i+1]-2*v[i]+v[i-1])/h^2 <= jU;</div>
<div>    </div>
<div># Boundary conditions for position</div>
<div>subject to x0:</div>
<div>    x[0] = 0;</div>
<div>subject to xf:</div>
<div>    x[N] = L;</div>
<div>    </div>
<div># Boundary conditions for velocity</div>
<div>subject to v0:</div>
<div>    v[0] = 0;</div>
<div>subject to vf:</div>
<div>    v[N] = 0;</div>
<div><br>
</div>
<div># Speed limit</div>
<div>subject to limit {i in 1..N-1}:</div>
<div>    v[i] <= if (x[i] <= 250) then 45/3.6 else if (x[i] <= 700) then 70/3.6 else 30/3.6;</div>
<div><br>
</div>
<div># Solve the optimization problem</div>
<div>solve;</div>
<div><br>
</div>
<div># write the data into a file for gnuplot</div>
<div>for {i in 0..N}</div>
<div>  printf : "%16.4e %16.4e %16.4e %16.4e \n", i*h, x[i], 3.6*v[i], 10*a[i] > ..\gnuplot.dat;</div>
<div><br>
</div>
<br>
<p></p>
</div>
</body>
</html>