2. A first example

In this section, a simple optimal control problem will be solved. Consider the optimal control problem for the Van der Pol oscillator model:

optimization VDP_Opt (objective = cost(finalTime),
                      startTime = 0,
                      finalTime = 20)

  // The states
  Real x1(start=0,fixed=true);
  Real x2(start=1,fixed=true);

  // The control signal
  input Real u;

  Real cost(start=0,fixed=true);

equation
  der(x1) = (1 - x2^2) * x1 - x2 + u;
  der(x2) = x1;
  der(cost) = x1^2 + x2^2 + u^2;
constraint 
   u<=0.75;
end VDP_Opt;

Create a new file named VDP_Opt.mop and save it in you working directory. Notice that this model contains both the dynamic system to be optimized and the optimization specification. This is possible since Optimica is an extension of Modelica and thereby supports also Modelica constructs such as variable declarations and equations. In most cases, however, Modelica models are stored separately from the Optimica specifications.

Next, create a Python script file and a write (or copy paste) the following commands:

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

# Import the plotting library
import matplotlib.pyplot as plt

Next, we transfer the model:

# Transfer the optimization problem to casadi
	op = transfer_to_casadi_interface("VDP_pack.VDP_Opt2", "VDP_Opt.mop")

The function transfer_to_casadi_interface invokes the CasADiInterface library to transfer the optimization problem and express it's components in terms of CasADi variables. This object represents the compiled model and is used to invoke the optimization algorithm:

res = op.optimize()

In this case, we use the default settings for the optimization algorithm. The result object can now be used to access the optimization result:

# Extract variable profiles
x1=res['x1']
x2=res['x2']
u=res['u']
t=res['time']

The variable trajectories are returned as numpy arrays and can be used for further analysis of the optimization result or for visualization:

plt.figure(1)
plt.clf()
plt.subplot(311)
plt.plot(t,x1)
plt.grid()
plt.ylabel('x1')

plt.subplot(312)
plt.plot(t,x2)
plt.grid()
plt.ylabel('x2')

plt.subplot(313)
plt.plot(t,u)
plt.grid()
plt.ylabel('u')
plt.xlabel('time')
plt.show()

You should now see the optimization result as shown in Figure 6.1, “Optimal profiles for the VDP oscillator”.

Figure 6.1. Optimal profiles for the VDP oscillator

Optimal profiles for the VDP oscillator

Optimal control and state profiles for the Van Der Pol optimal control problem.