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 compilation of models and the JMUModel class from pymodelica import compile_jmu from pyjmi import JMUModel # Import the plotting library import matplotlib.pyplot as plt
Next, we compile and load the model:
# Compile model jmu_name = compile_jmu("VDP_Opt","VDP_Opt.mop") # Load model vdp = JMUModel(jmu_name)
compile_jmu invokes the Optimica compiler and compiles the model into a DLL, which is then loaded when the
vdp object is created. This object represents the compiled model and is used to invoke the optimization algorithm:
res = vdp.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 8.1.