Permalink
Browse files

bump to casadi version 2.4.0

  • Loading branch information...
jgillis committed Aug 25, 2015
1 parent 14f677e commit 5f06bf9dff01d477a52653d3a54c2b7c5e5fce87
Showing with 192 additions and 145 deletions.
  1. +6 −7 examples/ocp.ipynb
  2. +2 −2 optoy/__init__.py
  3. +16 −3 optoy/base.py
  4. +46 −66 optoy/dynamic.py
  5. +19 −26 optoy/extensions/robustness.py
  6. +43 −0 optoy/simple_syntax.py
  7. +58 −39 optoy/static.py
  8. +2 −2 tests/test_dynamic.py
@@ -1,7 +1,7 @@
{
"metadata": {
"name": "",
"signature": "sha256:6125360fea9d1bac094f03709423a93ed03f9a484b0eb0237a60958098a77f7e"
"signature": "sha256:31e7b95b9eed04287708d52b0ed36129bfba3b318aa3dc9d4eaf254b408ff87a"
},
"nbformat": 3,
"nbformat_minor": 0,
@@ -25,21 +25,20 @@
"y.dot = x\n",
"q.dot = x**2+y**2\n",
"\n",
"ocp(T,[u>=-1,u<=1,q.start==0,x.start==1,y.start==0,x.end==0,y.end==0],T=T,N=20)"
"print ocp(T,[u>=-1,u<=1,q.start==0,x.start==1,y.start==0,x.end==0,y.end==0],T=T,N=20)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 16,
"output_type": "stream",
"stream": "stdout",
"text": [
"2.961575090066373"
"2.96157509007\n"
]
}
],
"prompt_number": 16
"prompt_number": 1
},
{
"cell_type": "code",
@@ -6,11 +6,11 @@
with open(os.path.join(os.path.dirname(__file__), 'conf.py')) as f:
exec(f.read())
from static import OptimizationVariable as var, OptimizationParameter as par, \
minimize, value
from dynamic import OptimizationState as state, OptimizationControl as control, ocp,\
time
from extensions.robustness import OptimizationDisturbance as dist, Prob, Sigma as cov
from simple_syntax import *
@@ -26,6 +26,17 @@
from casadi.tools import *
def try_expand(f):
""" Attempts to cast an MXFunction as SXFunxtion """
if not f.isInit():
f.init()
try:
r = SXFunction(f)
r.init()
return r
except:
return f
class OptimizationObject(MX):
"""
@@ -78,7 +89,7 @@ def cmpbyhash(self, b):
MX.__eq__ = cmpbyhash
# Get an exhausive list of all casadi symbols that make up f and gl
vars = set(getSymbols(veccat(el) if isinstance(el, list) else el))
vars = set(symvar(veccat(el) if isinstance(el, list) else el))
if dep:
while True:
@@ -167,8 +178,10 @@ def value(e, nums={}):
if e in OptimizationContext.eval_cache:
f, xp = OptimizationContext.eval_cache[e]
else:
syms = get_primitives(e, dep=False)
try:
syms = get_primitives(e, dep=False)
except:
return e
xp = []
for k in sorted(syms.keys()):
xp += syms[k]
@@ -46,22 +46,10 @@ def value_time(e, t):
try:
return DMatrix(e)
except:
f = MXFunction(getSymbols(e), [e])
f.init()
f.setInput(t)
f.evaluate()
return f.getOutput()
f = MXFunction("f",symvar(e), [e])
return f([t])[0]
def try_expand(f):
if not f.isInit():
f.init()
try:
r = SXFunction(f)
r.init()
return r
except:
return f
class OptimizationState(OptimizationContinousVariable):
@@ -109,9 +97,9 @@ def __init__(self, shape=1, lb=-inf, ub=inf, name="x", init=0):
def getDependent(cl, v):
newvars = set()
if hash(v) in cl._mapping:
newvars.update(set(getSymbols(cl._mapping[hash(v)].dot)))
newvars.update(set(symvar(cl._mapping[hash(v)].dot)))
if hash(v) in cl.lim_mapping:
newvars.update(set(getSymbols(cl.lim_mapping[hash(v)])))
newvars.update(set(symvar(cl.lim_mapping[hash(v)])))
return newvars
@@ -160,7 +148,7 @@ def __init__(self, shape=1, lb=-inf, ub=inf, name="u", init=0):
def getDependent(cl, v):
newvars = set()
if hash(v) in cl.lim_mapping:
newvars.update(set(getSymbols(cl.lim_mapping[hash(v)])))
newvars.update(set(symvar(cl.lim_mapping[hash(v)])))
return newvars
state = OptimizationState
@@ -227,6 +215,7 @@ def ocp(f, gl=[], regularize=[], verbose=False, N=20, T=1.0,
if not isinstance(gl, list):
raise Exception("Constraints must be given as a list")
f = f + OptimizationParameter()
print f.shape
# Determine nature of constraints, either g(x)<=0 or g(x)==0
gl_pure, gl_equality = sort_constraints(gl)
@@ -268,16 +257,14 @@ def ocp(f, gl=[], regularize=[], verbose=False, N=20, T=1.0,
for i in syms["p"]])
# Create the ode function
ode_out = MXFunction(syms["x"] +
ode_out = MXFunction("ode_out",syms["x"] +
syms["u"] +
syms["p"] +
syms["v"] +
ext_syms, [((T +
0.0) /
N) *
vertcat([i.dot for i in syms["x"]])])
ode_out.setOption("name", "ode_out")
ode_out.init()
# Group together all symbols that are constant over an integration interval
nonstates = struct_symMX([entry("u",
@@ -295,14 +282,13 @@ def ocp(f, gl=[], regularize=[], verbose=False, N=20, T=1.0,
ode_out_ins += nonstates[k, ...]
# Change the ode function signature to accept states/nonstates
ode = MXFunction(
daeIn(
x=states, p=nonstates), daeOut(
ode=ode_out(ode_out_ins)[0]))
ode.init()
ode = MXFunction("ode",
[states, nonstates], [
ode_out(ode_out_ins)[0]])
# Create the integrator
intg = explicitRK(ode, 1, 4, integration_intervals)
intg = simpleRK(ode, integration_intervals, 4)
intg = try_expand(intg)
Pw0 = P[...] + X["V", ...] + \
@@ -316,29 +302,26 @@ def ocp(f, gl=[], regularize=[], verbose=False, N=20, T=1.0,
ext_constr += ext_constr_e
# Set path constraints bounds
h_out = MXFunction(syms["x"] +
h_out = MXFunction("h_out",syms["x"] +
syms["u"] +
syms["p"] +
syms["v"] +
ext_syms, [a for a in gl_pure if dependsOn(a, syms["x"] +
syms["u"])])
h_out.setOption("name", "h_out")
g_out = MXFunction(syms["p"] +
ext_syms, [a for a in gl_pure if dependsOn(a, veccat(syms["x"] +
syms["u"]))])
g_out = MXFunction("g_out",syms["p"] +
syms["v"] +
ext_syms +
lims, [a for a in gl_pure if not dependsOn(a, syms["x"] +
syms["u"])])
g_out.setOption("name", "g_out")
f_out = MXFunction(syms["p"] + syms["v"] + ext_syms + lims, [f])
f_out.setOption("name", "f_out")
reg_out = MXFunction(syms["x"] +
lims, [a for a in gl_pure if not dependsOn(a, veccat(syms["x"] +
syms["u"]))])
f_out = MXFunction("f_out",syms["p"] + syms["v"] + ext_syms + lims, [f])
f_out.printDimensions()
reg_out = MXFunction("reg_out",syms["x"] +
syms["u"] +
syms["p"] +
syms["v"] +
ext_syms, [sumAll(vertcat([inner_prod(i, i) for i in regularize])) *
ext_syms, [sumRows(vertcat([inner_prod(i, i) for i in regularize])) *
T /
N])
reg_out.setOption("name", "reg_out")
# Expand if possible
h_out = try_expand(h_out)
@@ -347,7 +330,7 @@ def ocp(f, gl=[], regularize=[], verbose=False, N=20, T=1.0,
reg_out = try_expand(reg_out)
# Diagnostics
if dependsOn(f, syms["x"]):
if dependsOn(f, veccat(syms["x"])):
raise Exception(
"Objective function cannot contain pure state variables. Try adding .start or .end")
@@ -358,26 +341,25 @@ def ocp(f, gl=[], regularize=[], verbose=False, N=20, T=1.0,
[entry(str(i), expr=g) for i, g in enumerate(g_out(Pw0 + Lims))] +
[entry("path", expr=[h_out(X["X", k, ...] + X["U", k, ...] + Pw0) for k in range(N)])] +
ext_constr +
[entry("shooting", expr=[X["X", k + 1] - intg(x0=X["X", k], p=veccat([X["U", k]] + Pw0))[0] for k in range(N)])] +
[entry("shooting", expr=[X["X", k + 1] - intg([X["X", k], veccat([X["U", k]] + Pw0),1])[0] for k in range(N)])] +
([entry("periodic", expr=[X["X", -1] - X["X", 0]])]
if periodic else [])
)
# Build a regularization expression
# We dont use a helper state because we wisth to directly influence the
# We dont use a helper state because we wish to directly influence the
# objective
reg = sumAll(
reg = sumRows(
vertcat([reg_out(X["X", k, ...] + X["U", k, ...] + Pw0)[0] for k in range(N)]))
if reg.isempty(): reg = 0
# Construct the nlp
nlp = MXFunction(
nlp = MXFunction("nlp",
nlpIn(
x=X, p=P), nlpOut(
f=f_out(
Pw0 + Lims)[0] + reg, g=G))
nlp.setOption("name", "nlp")
nlp.init()
# Some extensions may wish to set default options
for cl in ext_cl:
exact_hessian = cl.setOptions(exact_hessian)
@@ -386,21 +368,19 @@ def ocp(f, gl=[], regularize=[], verbose=False, N=20, T=1.0,
if exact_hessian is None:
exact_hessian = True
# Allocate an ipopt solver
solver = NlpSolver("ipopt", nlp)
solver.setOption(
"hessian_approximation",
"exact" if exact_hessian else "limited-memory")
options = {}
options["hessian_approximation"] = "exact" if exact_hessian else "limited-memory"
if not verbose:
solver.setOption("print_time", False)
solver.setOption("print_level", 0)
solver.setOption("verbose", False)
solver.init()
options["print_time"] = False
options["print_level"] = 0
options["verbose"] = False
# Allocate an ipopt solver
solver = NlpSolver("solver", "ipopt", nlp, options)
# Set bounds on variables, set initial value
x0 = X(solver.input("x0"))
lbx = X(solver.input("lbx"))
ubx = X(solver.input("ubx"))
x0 = X(0)
lbx = X(0)
ubx = X(0)
for i in syms["v"]:
hs = str(hash(i))
@@ -420,19 +400,19 @@ def ocp(f, gl=[], regularize=[], verbose=False, N=20, T=1.0,
i.init, t=(k + 0.0) * T.init / N)
# Set parameter values
par = P(solver.input("p"))
par = P(0)
for i in syms["p"]:
h = str(hash(i))
par[h] = i.value
# Set constraint bounds
lbg = G(solver.input("lbg"))
ubg = G(solver.input("ubg"))
lbg = G(0)
ubg = G(0)
# Set normal constraints bounds
for i, eq in enumerate(
[e for g, e in zip(gl, gl_equality) if not dependsOn(g, syms["x"] + syms["u"])]):
[e for g, e in zip(gl, gl_equality) if not dependsOn(g,veccat(syms["x"] + syms["u"]))]):
if eq:
lbg[str(i)] = ubg[str(i)] = 0
else:
@@ -441,7 +421,7 @@ def ocp(f, gl=[], regularize=[], verbose=False, N=20, T=1.0,
# Set path constraints bounds
for i, eq in enumerate(
[e for g, e in zip(gl, gl_equality) if dependsOn(g, syms["x"] + syms["u"])]):
[e for g, e in zip(gl, gl_equality) if dependsOn(g, veccat(syms["x"] + syms["u"]))]):
if eq:
lbg["path", :, i] = ubg["path", :, i] = 0
else:
@@ -458,15 +438,15 @@ def ocp(f, gl=[], regularize=[], verbose=False, N=20, T=1.0,
lbg["periodic"] = ubg["periodic"] = 0
# Solve the problem numerically
solver.evaluate()
sol = solver(x0=x0,lbg=lbg,ubg=ubg,lbx=lbx,ubx=ubx,p=par)
# Raise an exception if not converged
if solver.getStat('return_status') != "Solve_Succeeded":
raise Exception(
"Problem failed to solve. Add verbose=True to see what happened.")
# Add the solution to the OptimizationObjects
opt = X(solver.output("x"))
opt = X(sol["x"])
# Extract solutions
for i in syms["v"]:
Oops, something went wrong.

0 comments on commit 5f06bf9

Please sign in to comment.