Annihilation of calcium oscillations (jModelica)

This page contains the source code to solve the Annihilation of calcium oscillations problem with JModelica. The automatic differentiation tool CasADI and the solver IPOPT were used to solve the problem. As this problem is numerically challenging, the following implementation does not result in a fully converged optimal solution.

Model file (calcium_opt.mop)

//-------------------------------------------------------------------------
//Annihilation of calcium oscillations with direct collocation using JModelica
//(c) Madeleine Schroter
//--------------------------------------------------------------------------


  optimization calcium_opt (objective = cost(finalTime),
                         startTime = 0,
                         finalTime = 22)
 
    // Parameters
    parameter Real k1 = 0.09;             // Parameter k1
    parameter Real k2 = 2.30066;          // Parameter k2
    parameter Real k3 = 0.64;             // Parameter k3
    parameter Real K4 = 0.19;             // Parameter K4
    parameter Real k5 = 4.88;             // Parameter k5
    parameter Real K6 = 1.18;             // Parameter K6
    parameter Real k7 = 2.08;             // Parameter k7
    parameter Real k8 = 32.24;            // Parameter k8
    parameter Real K9 = 29.09;            // Parameter K9
    parameter Real k10 = 5.0;             // Parameter k10
    parameter Real K11 = 2.67;            // Parameter K11
    parameter Real k12 = 0.7;             // Parameter k12
    parameter Real k13 = 13.58;           // Parameter k13
    parameter Real k14 = 153.0;           // Parameter k14
    parameter Real K15 = 0.16;            // Parameter K15
    parameter Real k16 = 4.85;            // Parameter k16
    parameter Real K17 = 0.05;            // Parameter K17
    parameter Real p17 = 1.0;             // Parameter p17
    parameter Real p18 = 1.0;             // Parameter p18
    parameter Real umax = 1.3;     
 
    parameter Real p1 = 100;              // Parameter p1
    parameter Real x0tilde = 6.78677;     // Parameter x0tilde 
    parameter Real x1tilde = 22.65836;    // Parameter
    parameter Real x2tilde = 0.384306;    // Parameter x2tilde 
    parameter Real x3tilde = 0.28977;     // Parameter x3tilde 
 
    // The states
    Real x0(start=0.03966, fixed=true);
    Real x1(start=1.09799, fixed=true);
    Real x2(start=0.00142, fixed=true);
    Real x3(start=1.65431, fixed=true);
    Real cost(start=0, fixed=true);
 
    // The control signal
    input Real u(free=true);
 
  equation
    der(x0) = k1+k2*x0-((k3*x0*x1)/(x0+K4))-((k5*x0*x2)/(x0+K6));
    der(x1) = k7*x0 - ((k8*x1)/(x1+K9));
    der(x2)=(k10*x1*x2*x3)/(x3+K11)+p18*k12*x1+k13*x0-(k14*x2)/(  p17*(1.0 + u*( umax - 1.0))  *x2+K15)-(k16*x2)/(x2+K17)+x3/10;
    der(x3)=-(k10*x1*x2*x3)/(x3+K11)+(k16*x2)/(x2+K17)-x3/10;
    der(cost)=(x0-x0tilde)*(x0-x0tilde)+(x1-x1tilde)*(x1-x1tilde)+
              (x2-x2tilde)*(x2-x2tilde)+(x3-x3tilde)*(x3-x3tilde)+p1*u;
 
  constraint
     x0>=0;
     x1>=0;
     x2>=0;
     x3>=0;

  end calcium_opt;

Run file

#Import the function for transfering a model to CasADiInterface
from pyjmi import transfer_optimization_problem

op=transfer_optimization_problem("calcium_opt", "calcium_opt.mop")

res=op.optimize()