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.

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.**

Option | Default | Description |
---|---|---|

`n_e` | 1 | Number of elements of the finite element mesh. |

`n_cp` | 20 | Number 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_options | None (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_phases | False | Specifies if the location of the phases should be allowed to be changed by the optimizer. Default: False |

`init_traj` | None (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_name` | Empty 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_phase`

s, `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.”.

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.

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

models.

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.”.

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.

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.