FAQ: how to add non CasADi code into a CasADi optimization?

See also: ⏯ CasADi callbacks, ⏯ External, 📖 User-defined function objects

The problem

For example, consider the case where you want to use gamma in your problem:

opti = casadi.Opti();

x = opti.variable();
opti.set_initial(x, 2);

opti.minimize(gamma(x));
opti.subject_to(gamma(sin(x))>=0.5);

opti.solver('ipopt');
sol = opti.solve();

(In Python, gamma is available with from math import gamma)

Now, this code will fail with RuntimeError: .../casadi/core/mx_node.cpp:385: 'to_double' not defined for class N6casadi10SymbolicMXE or Undefined function 'gamma' for input arguments of type 'casadi.MX'.

Why does it fail? Because CasADi developers have not implemented gamma to work with CasADi symbolic types. gamma only works with numbers, so you are stuck here.

The solution

The trick is to embed custom numerical-only code into a CasADi expression graph. You can do this by subclassing from casadi.Callback:

In Matlab, create a new file MyGamma.m:

classdef MyGamma < casadi.Callback
  methods
    function self = MyGamma()
      self@casadi.Callback();
      construct(self, 'Gamma', struct('enable_fd',true));
    end
    function [results] = eval(self, arg)
      % A cell array of DMs comes in
      arg = full(arg{1});

      % Some code that is not supported by CasADi
      res = gamma(arg);
      
      % A cell array of DMs should go out
      results = {res};
    end
  end
end

In Python, simply write the class-definition:

class MyGamma(casadi.Callback):
  def __init__(self):
    Callback.__init__(self)
    self.construct("Gamma", {"enable_fd": True})

  def eval(self, args):
    # A list of DMs comes in
    arg = np.array(args[0])

    # Some purely numeric code that is not supported by CasADi
    # Could also use if/while etc
    res = gamma(arg)

    # A list of DMs should go out
    results = [res]
    return results

Now instantiate this class: mygamma = MyGamma()

mygamma is now a CaADi Function mapping from a single scalar input to a single scalar output: Gamma:(i0)->(o0) CallbackInternal

You may proceed now to state your optimization as:

opti = casadi.Opti();

x = opti.variable();
opti.set_initial(x, 2);

opti.minimize(mygamma(x));
opti.subject_to(mygamma(sin(x))>=0.5);

opti.solver('ipopt');
sol = opti.solve();

This should solve.

Extras

More inputs

If your custom code needs more than a single input or output, you should add the method get_n_in/get_n_out to your class.

For example, if you would want 3 inputs, add the following to the class definition in Matlab:

function out = get_n_in(self)
   out = 3;
end

or Python:

def get_n_in(self):
   return 3

Matrix inputs

If your custom code needs matrix-valued inputs/outputs as opposed to scalar, you should add the method get_sparsity_in/get_sparsity_out

For example, if you want to work with a 3-by-3 matrix as first input, add the following to the class definition in Matlab:

function out = get_sparsity_in(self, i)
  if i==0 % Note: index-0 based
    out = casadi.Sparsity.dense(3,3);
  end
end

or Python:

def get_sparsity_in(self, i):
  if i==0:
    return casadi.Sparsity.dense(3,3)

Derivatives

We used the option enable_fd here to tell CasADi that it should use finite-differentiation on the mygamma function. For the other components of the CasADi expression graph (e.g. the sin(x)), CasADi will still be doing automatic differentiation. By overloading other Callback methods, you could also add exact derivatives for your piece of code.

Issues

You may get the error Callback object has been deleted if mygamma workspace variable goes out of scope or is overwritten before you solve the optimization. To avoid this, always keep a reference to mygamma (e.g. assign it to a variable with a unique name, make it a global variable, ...)