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 (objectiveIntegrand = x1^2 + x2^2 + u^2, startTime = 0, finalTime = 20) // The states Real x1(start=0,fixed=true); Real x2(start=1,fixed=true); // The control signal input Real u; equation der(x1) = (1 - x2^2) * x1 - x2 + u; der(x2) = x1; 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_optimization_problem # Import the plotting library import matplotlib.pyplot as plt
Next, we transfer the model:
# Transfer the optimization problem to casadi op = transfer_optimization_problem("VDP_pack.VDP_Opt2", "VDP_Opt.mop")
transfer_optimization_problem transfers the optimization problem into Python and expresses it's variables, equations, etc., using the automatic differentiation tool CasADi. 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”.