if_else usage

[Andres Codas Duarte]

Hi,

   I'm solving optimal control problems with ode's described by a set of equations which change depending on the state.  To write these equations in CasADi I use if_else.  This functionality works fine if the set of equations returns a valid number on all the domain.  But if I declare something like if_else(x>1,sqrt(x),x)  and if x goes negative the optimization solver brakes.  At the end of the message I provide you with a piece of code reproducing the problem.  Do you have any suggestion to overcome this?


   I'm also working with JModelica and I want to import statements involving if_else.  What is the syntax for this purpose?.  I tried in several ways (for example,  y = if x > 0 then x else -x; ) but I get always a compilation error when running compile_fmux.

   Thanks for your help!
  
   Andres

----------------------------------------------------------------------
from casadi import *

tf = 10.0
N = 1

ref = -2

x = ssym('x',2)
u = ssym('u',1)


dx0 = (x[1]-ref)*(x[1]-ref) + u*u
dx1 = if_else(x[1]>1,-sqrt(x[1]),-x[1])+u

fx = vertcat([dx0,dx1])

ode=SXFunction(daeIn(x=x,z=[],p=u,t=[]),daeOut(ode=fx))   

I = CVodesIntegrator(ode)
I.setOption("tf",tf/N)
I.setOption("verbose",True)
I.setOption('monitor',[ 'inputs' ,'outputs', 'res', 'resB', 'resQB' ,'reset', 'psetupB', 'djacB'])
I.init()                   

U = msym("U", N)

X0 = [0,2]

X  = msym(X0)
for k in range(N):
  (X,_,_,_) = I.call( (X,U[k]) )  

nlp = MXFunction(nlpIn(x=U),nlpOut(f=X[0],g=[]))

solver = IpoptSolver(nlp)
solver.init()


solver.setInput( 0.0,  "x0")
solver.solve()

u_opt = solver.output()

[joris.g...@gmail.com]

Dear Andres,

This is still an open issue in CasADi at the moment (issues https://github.com/casadi/casadi/issues/555 , https://github.com/casadi/casadi/issues/728):
In SX graphs, there is no proper branching support.

A quick fix for you is simply: if_else(x[1]>1,-sqrt(fabs(x[1])),-x[1])+u

Best regards,
   Joris

[Andres Codas Duarte]

Dear Joris,

   Thanks for your answer.  Your suggestion works fine for my example.  But still I need some help to implement it.  As I said, I'm importing the model with modelica.

   It seems that there is no way to import 'fabs' from Modelica, since the parsing function "SX CasADi::SymbolicOCP::readExpr" does not implement it.

   A workaround would be to pick all expressions involving sqrt(...) and replace them with (sqrt(fabs(...))).  How can I make this substitutions?

   There is also no way to import "if_else" but it seems you have a workaround using "NoEvent", but I couldn't figure out the syntax for make it work...  How do I write NoEvent(cond,exp_True,exp_False) in modelica to be able to parse it? 

    A side comment:  I couldn't find how you implement fabs, but this is in a way 'branching' since you have to do if_else(x>0,x,-x) ... so, it is a bit surprising for me the workaround you proposed.

   Thanks again!

   Andres

[Joel Andersson]

A comment on this: fabs(x) is a non-smooth function but it does not involve any branching. It has a derivative rule which is "almost everywhere" equal to sign(x). sign(x) in turn has a derivative which is almost everywhere equal to zero. Therefore it works just fine in an algorithmic differentiation context, if you keep the "almost everywhere" in mind.

About if_then_else, the corresponding atomic operation is something like "if_else_zero", defined like: f(x,y) := x ? y : 0. This function is easier to handle in an algorithmic differentiation context since it has a simple derivative rule df/fx = 0 and df/dy = x, again in the "almost everywhere" sense.

CasADi then defines if_then_else as:

f(cond,x,y) := if_else_zero(cond, x) + if_else_zero(not cond, y)

I think it is defined in sx_tools.hpp. It should be included when you include symbolic/casadi.hpp.

I hope this made things clearer!

Best,

Joel