7. Optimization of ODEs using Pseduo-Spectral methods

The automatic differentiation tool CasADi is connected to IPOPT. The collocation algorithms are the Gauss Pseudospectral-, the Legendre Pseudospectral- and the Radau Pseudospectral-method. These methods can be used either as a local method or a global method where in the local case the time horizon is discretized into elements and each element is approximated with a polynomial. In the global case, a polynomial is used over the complete time horizon. For a mathematic overview of the methods, see the Python API documentation for the class PseudoSpectral.

This section relies on a model format called FMUX, and the Python class CasadiModel. These are supported in JModelica.org, but no longer used for the main CasADi-based optimization tool chain. For more information about the FMUX format and the CasadiModel class, see the JModelica.org-1.12 User Guide, available from http://jmodelica.org/page/236.

The implementation resides in the module casadi_interface and contains a model class CasadiModel. A compiler method for creating models compliant with CasadiModel resides in pymodelica, compile_fmux. Usage of the methods is consistent with the work flow for optimization of JMUModels. Shown below is a Python code example describing how to invoke the optimization method.

from pymodelica import compile_fmux
from pyjmi import CasadiModel

casadi_name = compile_fmux("MyModel", "MyModels.mop")
model = CasadiModel(casadi_name)

res = model.optimize(algorithm="CasadiPseudoSpectral")

#Plot...

As can be seen from the above code example the real difference from the native optimization method is that here we use a different model class and a different compilation class, except from that, it can be used in similar fashion.

7.1. Options

Options can be viewed and set via the options dictionary which can then be passed to the optimize call. The option dictionary can be retrieved with the below code example where model is a CasadiModel instance. It is also shown how to use the interactive help.

opts = model.optimize_options(algorithm="CasadiPseudoSpectralAlg")

opts? #View the help text

To provide the options to the optimize call, use the options attribute.

res = model.optimize(algorithm="CasadiPseudoSpectralAlg", options=opts)

In Table 7.4, “Options for the Pseudospectral optimization algorithms.” the options for the PseudoSpectral methods are shown.

Table 7.4. Options for the Pseudospectral optimization algorithms.

OptionDefaultDescription
n_e1Number of elements of the finite element mesh.
n_cp20Number of collocation points in each element. Values between 1 and 80 are supported.
discr"LG"Determines the discretization of the problem. Can be either "LG" (Legendre-Gauss), LGR (Legendre-Gauss-Radau), LGL (Legendre-Gauss-Lobatto).
link_options[] (not used)This option allows users to specify states that are allowed to be discontinuous between phases (elements) and to connect the transition with a model parameter. Example, [(1,"x1","dx1")], this allowes the variable x1 to be discontinuous between phase 1 and 2 with the parameter dx1. The generating constraint becomes, x1_N^1 - x1_0^2 - dx1 = 0. There is no limit that the same parameter can be used in multiple transition such as, [(1,"x1","dx1"), (1,"x2","dx1"), (2,"x1","dx1")]
phase_optionsNone (not used)This options allows users to connect parameters to phase boundaries (time). Example, in a three phase problem the parameters t1 and t2 can be specified to be the boundaries of the phases such as, ["t1", "t2"], the option free_phases have also be set to true. Default: None
free_phasesFalseSpecifies if the location of the phases should be allowed to be changed by the optimizer. Default: False
init_trajNone (i.e. not used, set this argument to activate initialization)Variable trajectory data used for initialization of the optimization problem. The data is represented by an object of the type pyjmi.common.io.ResultDymolaTextual.
result_file_nameEmpty string (default generated file name will be used)Specifies the name of the file where the optimization result is written. Setting this option to an empty string results in a default file name that is based on the name of the optimization class.
result_format'txt'Specifies in which format to write the result. Currently only textual mode is supported.


Options that need a little more attention are the free_phases, link_options and phase_options. These options allows a user to setup a problem involving discontinuities. As an example, say that in a model a state is allowed to be changed discontinuous at an arbitrary time point, tp, with the parameter dv. See Figure 7.14, “Handling of discontinuities.”.

Figure 7.14. Handling of discontinuities.

Handling of discontinuities.

In this case we set up the solver to use two elements (phases) and at the phase boundary we impose an extra constraint that our state in phase one plus the discontinuity should be equal to our state in phase two.

opts = model.optimize_options(algorithm="CasadiPseudoSpectralAlg")

opts["n_e"]=2 #Two elements
opts["link_options"] = [(1,"v","dv")]   # The parameter dv should be connected
                                        # between phase 1 and 2 to variable v
opts["phase_options"] = ["t1"]          # The parameter t1 connects phase 1 and 2
opts["free_phases"] = True              # The phase boundary is allowed to be changed.

Note that the parameters dv and t1 should be defined in the model and should be set to free.

7.2. Examples

In the following subsections it will be shown how to use the PseudoSpectral methods on Optimica models.

7.2.1. Van der Pol oscillator

This example shows the well known Van der Pol (VDP) oscillator being optimized as an optimal control problem. The formulation is somewhat different than in previous example of the VDP problem, as here there is no extra cost state, instead the cost function have been directly defined in the Lagrange term (objectiveIntegrand). There is also no bounds on u. Below, the Optimica for the example is shown.

optimization VDP_Opt (objectiveIntegrand = (x1^2 + x2^2 + u^2),
                         startTime = 0,
                         finalTime = 20)

    // State start values
    parameter Real x1_0 = 0;
    parameter Real x2_0 = 1;

    // The states
    Real x1(start = x1_0, fixed=true);
    Real x2(start = x2_0, fixed=true);

    // The control signal
    input Real u;

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

To compile and load the VDP model into JModelica.org we first import the necessary Python modules. Then we call the compiler method and we load it into our model class. Finally the optimization is performed by a call to the optimize method. An extra note in this example is that we change the number of collocation points from the default 20 to 40 via the options dictionary.

from pymodelica import compile_fmux
from pyjmi import CasadiModel

casadi_name = compile_fmux("VDP_Opt", "VDP.mop")
model = CasadiModel(casadi_name)

# Get the options dictionary
opts = model.optimize_options(algorithm="CasadiPseudoSpectralAlg")
opts["n_cp"] = 40 # Number of collocation points

res = model.optimize(algorithm="CasadiPseudoSpectralAlg", options=opts)

#Plot...

The result can then be used to visualize the result and a resulting figure is shown in Figure 7.15, “The Van der Pol oscillator optimized using Gauss PseudoSpectral method.”.

Figure 7.15. The Van der Pol oscillator optimized using Gauss PseudoSpectral method.

The Van der Pol oscillator optimized using Gauss PseudoSpectral method.

7.2.2. Hohmann Transfer

This example will show how to use the pseudo-spectral method for solving a problem of altitude change of a satellite orbiting earth. The problem have a well known solution called a Hohmann Transfer from the German scientist Walter Hohmann. The solution consists of two rocket burns, one for leaving the current orbit and one when the final orbit have been reached, both perpendicular to the radial component in a rotation coordinating system. The magnitude of the rocket burns can also be calculated theoretically.

We will solve this problem using the gauss pseudo-spectral method with two phases. The phases will be connected in such a way that the velocities are allowed to be discontinuous, see Figure 7.14, “Handling of discontinuities.” . The magnitude of the discontinuity will correspond to the first rocket burn. The second rocket burn are connected as constraints at the final time.

We start by setting up our model of a planar two-body system where the bodies are the earth and the satellite. We also scale the problem with the orbital time and the radius of the earth.

model Orbital
  //"State variables"
  Real x(fixed=true,start=x0);
  Real y(fixed=true,start=y0);
  Real vx(fixed=true,start=vx0);
  Real vy(fixed=true,start=vy0);
  
  //"Initial values"
  parameter Real x0=rStart;
  parameter Real y0=0;
  parameter Real vx0=0;
  parameter Real vy0=vStart;
  
  //"Constants"
  constant Real G=6.67384e-11 "Gravitational Constant";
  constant Real M=5.9742e24 "Earth Mass";
  constant Real R=6.371e6 "Earth Radii";
  
  //"Parameters"
  parameter Real my=G*M/R^3*T^2;
  parameter Real T=5400; //Orbital time of rStart

  parameter Real rFinal = rStart*2 "Radii of final orbit";
  parameter Real rStart = (R+300e3)/R "Radii of start orbit";
  parameter Real vStart = 7730/R*T "Velocity of start orbit";
  
equation //Equation of motion of a planar two-body problem
  der(x)=vx;
  der(y)=vy;
  der(vx)=-x*my/(x^2+y^2)^(3/2);
  der(vy)=-y*my/(x^2+y^2)^(3/2);
end Orbital;

Our starting orbit is 300 km above the earths surface and our target orbit is 6971 km.

Next, we set up our optimization problem.

optimization HohmannTransfer(objective=(dx1^2+dy1^2)^0.5+(dx2^2+dy2^2)^0.5, startTime=0, finalTime(free=true,min=0.5,max=10,initialGuess=0.8))
  extends Orbital;

  //"Rocket burns"
  parameter Real dx1(free=true,min=-2,max=2,start=0.05);
  parameter Real dy1(free=true,min=-2,max=2,start=0.05);
  parameter Real dx2(free=true,min=-2,max=2,start=0.05);
  parameter Real dy2(free=true,min=-2,max=2,start=0.05);
  
  //"Time for firing the first rockets"
  parameter Real t1(free=true,min=0.001,max=0.01, start=0.005);

constraint
  (x(finalTime)^2+y(finalTime)^2)^0.5=rFinal; //"Final position constraint"
  (vx(finalTime)+dx2)^2+(vy(finalTime)+dy2)^2=my/rFinal; //"Final velocity"
  x(finalTime)*(vx(finalTime)+dx2)+y(finalTime)*(vy(finalTime)+dy2)=0; //"Final Radial velocity component"
end HohmannTransfer;

As can be seen from the above the code, our objective is to minimize the rocket burns (d1 and d2) with the final time unknown. The constraints are set on the final position (circular) and on the final velocities so that we are actually staying in the target orbit. There is also a parameter t1 set to be free which will be connected to the problem via Python and will be used to specify the time of the first rocket burn. From the problem, one can also see that dx1 and dy1 are not included in the problem, except for the objective. They will be connected via Python and determine the velocity changes between phase 1 and phase 2. Note that in the following scripts the models have been stored as Hohmann.mop .

We start our script by importing the necessary modules and methods for compiling the two different models and also for loading the compiled models into the JModelica.org framework. Note that it is necessary to disable the compiler option normalize_minimum_time_problems when solving problems with free start or final time using this algorithm.

import numpy as N
import matplotlib.pyplot as plt

from pymodelica import compile_fmu
from pyfmi import load_fmu

from pymodelica import compile_fmux
from pyjmi import CasadiModel

jmu_name  = compile_fmu("Orbital","Hohmann.mop")
comp_opts = {"normalize_minimum_time_problems": False}
fmux_name = compile_fmux("HohmannTransfer", "Hohmann.mop", compiler_options=comp_opts)

Next, the optimization model needs to be loaded and options specified.

# Optimization
model = CasadiModel(fmux_name)
    
opts = model.optimize_options(algorithm="CasadiPseudoSpectralAlg")

opts["n_cp"] = 40                                      # Number of collocation points
opts["n_e"] = 2                                        # Number of phases
opts["free_phases"] = True                             # The phase boundary is allowed to be changed in time
opts["phase_options"] = ["t1"]                         # The phase boundary is connected to variable t1
opts["link_options"] = [(1,"vy","dy1"),(1,"vx","dx1")] # Allow for discontinuities between phase 1 and 2 for vy and vx.
                                                       # The discontinuities are connected by dy1 and dx1

In the above code we have loaded our optimization model fmux_name into a class, CasadiModel. Our options dictionary is retrieved by calling the optimization options method with our method of choice as input, CasadiPseudoSpectral. We then specify the number of collocation points in each phase and also the number of phases. The next options are a bit more interesting. The option free_phases specifies that the boundary between phase 1 and phase 2 are allowed to be changed by the optimizer and the phase_options specifies that the boundary should be connected to the parameter t1 in the model where we also have specified constraints on the parameter. Finally we connect the first rocket burn (dx1 and dy1) to the corresponding variables, dx and dy, allowing them to be changed on the phase boundary.

Now we are ready to optimize our problem.

# Optimize
res_opt = model.optimize(algorithm="CasadiPseudoSpectralAlg",options=opts)
    
# Get results
dx1,dy1,dx2,dy2 = res_opt.final("dx1"), res_opt.final("dy1"), res_opt.final("dx2"), res_opt.final("dy2")
r1,r2,my        = res_opt.final("rStart"), res_opt.final("rFinal"), res_opt.final("my")
tf,t1           = res_opt.final("time"), res_opt.final("t1")

The above code performs the optimization and also retrieves the values used to verify the solution.

To verify the solution we perform simulations using the resulting rocket burns and even simulate further in time to verify that we are staying in the targeted orbit.

# Verify solution by simulation
model = load_fmu(fmu_name)
    
# Simulation of Phase 1
res = model.simulate(final_time=t1,options={'ncp':100,'solver':"CVode"})

x_phase_1,y_phase_1 = res["x"], res["y"]
    
# Simulation of Phase 2
model.set("vx",dx1+res["vx"][-1])
model.set("vy",dy1+res["vy"][-1])

res = model.simulate(start_time=t1,final_time=tf,options={'ncp':100,'solver':"CVode","initialize":False})

x_phase_2,y_phase_2 = res["x"], res["y"]

# Simulation of Phase 3 (verify that we are still in orbit)
model.set("vx",dx2+res["vx"][-1])
model.set("vy",dy2+res["vy"][-1])

res = model.simulate(start_time=tf,final_time=tf*2,options={'ncp':100,'solver':"CVode","initialize":False})

x_phase_3,y_phase_3 = res["x"], res["y"]

In the above simulation code we used the model compiled in the beginning of the script ,"Orbital" and then performed three simulations. One for the first phase, one for the second phase and finally an extra simulation for verification. In between the simulations we added the calculated rocket burns to the velocities and also retrieved the result.

Next, we set up the figure where we plot the the initial and target orbit together with the earth and the satellite trajectory. We also plot the rocket burns.

# Plot Earth
r = 1.0
xE = r*N.cos(N.linspace(0,2*N.pi,200))
yE = r*N.sin(N.linspace(0,2*N.pi,200))
plt.plot(xE,yE,label="Earth")

# Plot Orbits
r = res.final("rStart")
xS = r*N.cos(N.linspace(0,2*N.pi,200))
yS = r*N.sin(N.linspace(0,2*N.pi,200))
plt.plot(xS,yS,label="Low Orbit")

r = res.final("rFinal")
xSE = r*N.cos(N.linspace(0,2*N.pi,200))
ySE = r*N.sin(N.linspace(0,2*N.pi,200))
plt.plot(xSE,ySE,label="High Orbit")
        
# Plot Satellite trajectory
x_sim=N.hstack((N.hstack((x_phase_1,x_phase_2)),x_phase_3))
y_sim=N.hstack((N.hstack((y_phase_1,y_phase_2)),y_phase_3))

plt.plot(x_sim,y_sim,"-",label="Satellite")
        
# Plot Rocket Burns
plt.arrow(x_phase_1[-1],y_phase_1[-1],0.5*dx1,0.5*dy1, width=0.01,label="dv1")
plt.arrow(x_phase_2[-1],y_phase_2[-1],0.5*dx2,0.5*dy2, width=0.01,label="dv2")
        
plt.legend()
plt.show()

In Figure 7.16, “Hohmann Transfer” the result is visualized. The arrows represents the rocket burns.

Figure 7.16. Hohmann Transfer

Hohmann Transfer

7.3. Limitations

Current limitations for the PseudoSpectral methods include,

  • Final time or start time cannot be explicitly in the cost function. Although states of final time /start time can be used.

  • No path constraints. However variable bounds can be used.

  • Only explicit ordinary differential equations.

  • Still under development and not extensively tested, consider using JMUModel.optimize as a first choice.

  • No support for scaling.

  • Parameters cannot be set after compilation.