Setup of Algebraic equations in DAE

bransfor
2012-05-16
2013-04-11
  • bransfor
    bransfor
    2012-05-16

    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()
    
     
  • Joel Andersson
    Joel Andersson
    2012-05-21

    Hello Bransfor,

    In the current trunk version of CasADi, users are expected to provide state derivatives of algebraic states, though of course they will never enter in the DAE residual equations. The way algebraic states are handled in CasADi right now is not satisfactory, and it creates problems in particular for the adjoint sensitivity analysis.

    There is a branch of CasADi, branches/new_integrator2, where this has been changed and the differential states have been separated from the algebraic states. Also options like "is_differential" disappear. This branch is working and in sync with the trunk, and will eventually be merged with the trunk, and when this happens, there are some syntax changes for the users. The only reason that it hasn't happened yet, is that the work on the adjoint sensitivity analysis isn't finished. Anyway, it should already work just as well (or better) than the trunk.

    As for the options, they do not propagate from the Simulator class to the integrator class(es). When the re factoring of the integrator class has been completed, I'll take a look at the Simulator class. It is possible that its functionality can be merged into the integrator base class, though I'm not sure about that yet.

    Best,
    Joel

     
  • Joel Andersson
    Joel Andersson
    2012-10-08

    As mentioned, algebraic variables are now separated from differential states. Algebraic variables are now an entirely internal issue for the integrators.