Trying migration for AMPL to Pyomo

[pablope...@gmail.com]

Hi everyone,

I came from ampl world, but recently I have been triying to migrate to open source technologies. After testing Casadi, i understand that is weird syntax and its learning curve it was not for me, so I give a trial to Pyomo. I found it great, but maybe because I am kind of new in Python language and/or because I am su used to the ampl path, I am not even able of doing something more complex that the simple examples adjoined to Pyomo.

The test I have been trying to do is migrating the ampl car sample that cames with Ipopt source code distribution into a Pyomo file, but i am not able of make it working:

I have problems with the end condition (reach a location with zero speed at final time) and with the cost function (minimize final time, where time is the continuos set). Please, might someone try helping me with this two issues?

Thanks in advance. Here is the code:

# min tf

# dx/dt = 0

# dv/dt = a - R*v^2

# x(0) = 0; x(tf) = 100

# v(0) = 0; v(tf) = 0

# -3 <= a <= 1 (a is the control variable)

#!Python3.5

from pyomo.environ import *

from pyomo.dae import *

N = 20;

T = 10;

L = 100;

m = ConcreteModel()

# Parameters

m.R = Param(initialize=0.001)

# Variables

def x_init(m, i):

    return i*L/N

m.t = ContinuousSet(bounds=(0,1000))

m.x = Var(m.t, bounds=(0,None), initialize=x_init)

m.v = Var(m.t, bounds=(0,None), initialize=L/T)

m.a = Var(m.t, bounds=(-3.0,1.0), initialize=0)

# Derivatives

m.dxdt = DerivativeVar(m.x, wrt=m.t)

m.dvdt = DerivativeVar(m.v, wrt=m.t)

# Objetives

m.obj = Objective(expr=m.t[N])

# DAE

def _ode1(m, i):

    if i==0:

        return Constraint.Skip

    return m.dxdt[i] == m.v[i]

m.ode1 = Constraint(m.t, rule=_ode1)

def _ode2(m, i):

    if i==0:

        return Constraint.Skip

    return m.dvdt[i] == m.a[i] - m.R*m.v[i]**2

m.ode2 = Constraint(m.t, rule=_ode2)

# Constraints

def _init(m):

    yield m.x[0] == 0

    yield m.v[0] == 0

    yield ConstraintList.End

m.init = ConstraintList(rule=_init)

'''

def _end(m, i):

    if i==N:

        return m.x[i] == L amd m.v[i] == 0

    return Constraint.Skip

m.end = ConstraintList(rule=_end)

'''

# Discretize

discretizer = TransformationFactory('dae.finite_difference')

discretizer.apply_to(m, nfe=N, wrt=m.t, scheme='BACKWARD')

# Solve

solver = SolverFactory('ipopt', executable='C:\\EXTERNOS\\COIN-OR\\win32-msvc12\\bin\\ipopt')

results = solver.solve(m, tee=True)

[Bethany Nicholson]

I've implemented your example using pyomo.dae and added it to the examples in the github repository. Please see here. The key to implementing this problem is to reformulate it by removing the time scaling from the ContinuousSet. In the example notice the difference between m.tau and m.time. Also notice that the differential equations are scaled by the final time value. Let me know if you have any questions about this reformulation.

Bethany

--
You received this message because you are subscribed to the Google Groups "Pyomo Forum" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pyomo-forum+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Santiago rodriguez]

Great example.

Santiago

To unsubscribe from this group and stop receiving emails from it, send an email to pyomo-forum...@googlegroups.com.

[pablope...@gmail.com]

Thank you for your soon answer and the nice reimplementation of the sample.

As I have told, i come from ampl, so i am used to automatic differentiation but not so used to automatic discretization. Also Pyomo systax is nice, but a little more verbose than the ampl one.

About your solution, I still have a few questions that you may help me to understand:

1) Is it necessary to protect ODEs and constraints with the "if i==0 return constraint.skip" clause? I am not looking any problem in use a general case for all the elements, but probably there must be something I have not had into account.

2) As everybody who has worked with nlp and optimal control knows, interior point methods ar so sensitive to the initial conditions set in the variables, being sometimes this the difference between an infeasible problem and an optimal solution found. Which is the right way of initialize the variables, per example "var x{i in 0..N} >= 0, := i*(L/N);" using an ampl syntax?

3) Related with #2, ampl allows warm start, where you can use the results of a system resolution as initial values for a new optimization problem. Does Pyomo allow this feature in an easy way?

Thanks a lot,

Pablo

[Bethany Nicholson]

1) Those lines are necessary right now to keep things completely consistent. However, you'll likely find that the problem still solves if you remove those lines. We're working on an update to the DAE package which will do that constraint skipping for you depending on the discretization you select.

2) If you wanted to initialize x like that it would be:

def _initX(m,t):

    return t*m.L

m.x = Var(m.tau, bounds=(0,None), initialize=_initX)

(I left out the N because I'm assuming tau is bounded between 0 and 1)

3) You can also warm start the problem in Pyomo. After you solve the problem the solution is loaded back into the model so if you solve it again in the same script your model will be warm started. You can also save the results from a previous run in a Python dictionary and manually load that data and set the values on your variables. If you want to save the multipliers you can get those using suffixes similar to ampl. To summarize, there are ways to warm start your model but it is up to you to implement them (figure out what values to save, how to store them, and how to load them back into your model)

Bethany

--
You received this message because you are subscribed to the Google Groups "Pyomo Forum" group.

To unsubscribe from this group and stop receiving emails from it, send an email to pyomo-forum+unsubscribe@googlegroups.com.