My questions is prompted by this topic. Does one need to specify a differential for every state when creating a DAE instance? For example, I've coded example 9.3 in L. Biegler's nonlinear programming book. The third equation is algebraic, so I did not include a differential for 'yC'.

Also, does the options set in the Integrator class propagate to the simulator class?

The code works and produces a figure resembling Figure 9.3 in the book. However I would still like to know if there is anything wrong with this approach that might cause problem when I try to use it for other DAEs.

import casadi import numpy import pylab calc_ic = 0 # independent variable time = casadi.SX('time') # decision variables zA = casadi.SX('zA') zB = casadi.SX('zB') yC = casadi.SX('yC') # differentials dzA = casadi.SX('dzA/dt') dzB = casadi.SX('dzB/dt') # parameters p1 = casadi.SX('p1') p2 = casadi.SX('p2') p3 = casadi.SX('p3') p4 = casadi.SX('p4') # residual error res = [(-p1 * zA + p2 * zB) - dzA, (p1 * zA - (p2 + p3) * zB + p4 * yC) - dzB, (zA + zB + yC) - 1] # construct DAE input dae_in = [[]] * casadi.DAE_NUM_IN dae_in[casadi.DAE_T] = [time] dae_in[casadi.DAE_Y] = [zA, zB, yC] dae_in[casadi.DAE_YDOT] = [dzA, dzB] dae_in[casadi.DAE_P] = [p1, p2, p3, p4] # construct function ffcn = casadi.SXFunction(dae_in, [res]) # instantiate integrator integrator = casadi.IdasIntegrator(ffcn) integrator.setOption('calc_ic', calc_ic) integrator.setOption('is_differential', [1, 1, 0]) integrator.init() # instantiate Simulator time = numpy.linspace(0, 1) simulator = casadi.Simulator(integrator, time) simulator.init() simulator.setInput([1, 0, 0], casadi.INTEGRATOR_X0) simulator.setInput([3.997, 1.998, 40.538, 20.264], casadi.INTEGRATOR_P) if not calc_ic: simulator.setInput([0, 0, 0], casadi.INTEGRATOR_XP0) # run simulator and collect trajectory simulator.evaluate() output = simulator.output().toArray() # plot trajectory ax = pylab.subplot(111) ax.plot(time, output) ax.set_xlabel('time') ax.set_ylabel('Normalized Concentration') pylab.legend(('zA', 'zB', 'yC')) pylab.show()