Next, open a Python shell, preferably using the Pylab mode. If you are running Windows, select the menu option provided with the JModelica.org installation. If you are running Linux or Mac OS X, open a terminal and enter the command:

> /path/to/jmodelica_installation/Python/jm_ipython.sh -pylab

As your first action, go to the working directory you have created:

In [1]: cd '/path/to/working/directory'

In order to run the Python script, use the 'run' command:

in [2]: run cstr.py

The functions and classes used in the tutorial script need to be imported into the Python script. This is done by the following Python commands. Copy them and past them either directly into you Python shell or, preferably, into your Python script file.

import os.path
from jmodelica import initialize
from jmodelica import simulate
from jmodelica import optimize

import jmodelica.jmi as jmi
from jmodelica.compiler import OptimicaCompiler

import numpy as N
import matplotlib.pyplot as plt

Before we can do operations on the model, such as optimizing it, the model file must be compiled and the resulting DLL file loaded in Python. These steps are described in more detail in the tutorial on Compilation of models.

# Create a Modelica compiler instance
oc = OptimicaCompiler()
    
# Compile the stationary initialization model into a DLL and load it
init_model = oc.compile_model("CSTR.CSTR_Init", "CSTR.mo", target='ipopt')

At this point, you may open the file CSTR.mo, containing the CSTR model and the static initialization model used in this section. Study the classes CSTR.CSTR and CSTR.CSTR_Init and make sure you understand the models. Before proceeding, have a look at the interactive help for one of the functions you used:

In [8]: help(oc.compile_model)

In the next step, we would like to specify the first operating point, A, by means of a constant input cooling temperature, and then solve the initialization problem assuming that all derivatives are zero.

# Set inputs for Stationary point A
Tc_0_A = 250
init_model.set_value('Tc',Tc_0_A)
    
# Solve the DAE initialization system with Ipopt
(init_model, init_result) = initialize(init_model)
    
# Store stationary point A
c_0_A = init_result.get_variable_data('c').x[0]
T_0_A = init_result.get_variable_data('T').x[0]

# Print some data for stationary point A 
print(' *** Stationary point A ***') 
print('input Tc = %f' % Tc_0_A) 
print('state c = %f' % c_0_A) 
print('state T = %f' % T_0_A)

Notice how the function set_value is used to set the value of the control input. The initialization algorithm is invoked by calling the function 'initialize', which returns the initialization result in the object 'init_result'. The 'initialize' function relies on the algorithm Ipopt for computing the solution of the initialization problem. The values of the states corresponding to grade A can then be extracted from the result object. Look carefully at the printouts in the Python shell to see a printout of the stationary values. Display the help text for the 'initialize' function and take a moment to look through it. The procedure is now repeated for operating point B:

# Set inputs for Stationary point B
Tc_0_B = 280
init_model.set_value('Tc',Tc_0_B)
        
# Solve the DAE initialization system with Ipopt
(init_model, init_result) = initialize(init_model)
    
# Store stationary point B
c_0_B = init_result.get_variable_data('c').x[0]
T_0_B = init_result.get_variable_data('T').x[0]

# Print some data for stationary point B 
print(' *** Stationary point B ***') 
print('input Tc = %f' % Tc_0_B) 
print('state c = %f' % c_0_B) 
print('state T = %f' % T_0_B)

We have now computed two stationary points for the system based on constant control inputs.

The optimal control problem we are about to solve is given by:

and is expressed in Optimica format in the class CSTR.CSTR_Opt in the CSTR.mo file above. Have a look at this class and make sure that you understand how the optimization problem is formulated and what the objective is.

Direct collocation methods often require good initial guesses in order to ensure robust convergence. Since initial guesses are needed for all discretized variables along the optimization interval, simulation provides a convenient mean to generate state and derivative profiles given an initial guess for the control input(s). It is then convenient to set up a dedicated model for computation of initial trajectories. In the model CSTR.CSTR_Init_Optimization in the CSTR.mo file, a step input is filtered through a first order filter in order to generate a smooth input for the CSTR system. The filtering is done in order not to excite unstable modes of the system, and in particular to avoid sudden ignition. Notice also that the variable names in the initialization model must match those in the optimal control model. Therefore, also the cost function is included in the initialization model.

Start by creating an input trajectory to be passed to the simulator:

# Create the time vector
t = N.linspace(1,150.,100)
# Create the input vector from the target input value. The
# target input value is here increased in order to get a
# better initial guess.
u = (Tc_0_B+35)*N.ones(N.size(t,0))
# Create a matrix where the first column is time and the second column represents
# the input trajectory.
u_traj = N.transpose(N.vstack((t,u)))

Next, compile the model and set model parameters:

# Compile the optimization initialization model and load the DLL
init_sim_model = oc.compile_model("CSTR.CSTR_Init_Optimization", "CSTR.mo", target='ipopt')

# Set model parameters
init_sim_model.set_value('cstr.c_init',c_0_A)
init_sim_model.set_value('cstr.T_init',T_0_A)
init_sim_model.set_value('Tc_0',Tc_0_A)
init_sim_model.set_value('c_ref',c_0_B)
init_sim_model.set_value('T_ref',T_0_B)
init_sim_model.set_value('Tc_ref',u[0])

Having initialized the model parameters, we can simulate the model using the 'simulate' function.

(init_sim_model,res) = simulate(init_sim_model,alg_args={'start_time':0.,'final_time':150.,
                                                         'input_trajectory':u_traj})

The function 'simulation' first computes consistent initial conditions and then simulates the model in the interval 0 to 150 seconds with the input trajectory specified by 'u_traj'. Notice that the arguments to the simulation function is specified in a Python dictionary. Take a moment to read the interactive help for the 'simulate' function.

The simulation result is returned in the output argument 'res', from which you may now retrieve trajectories for plotting:

# Extract variable profiles
c_init_sim=res.get_variable_data('cstr.c')
T_init_sim=res.get_variable_data('cstr.T')
Tc_init_sim=res.get_variable_data('cstr.Tc')

# Plot the results
plt.figure(1)
plt.clf()
plt.hold(True)
plt.subplot(311)
plt.plot(c_init_sim.t,c_init_sim.x)
plt.grid()
plt.ylabel('Concentration')

plt.subplot(312)
plt.plot(T_init_sim.t,T_init_sim.x)
plt.grid()
plt.ylabel('Temperature')

plt.subplot(313)
plt.plot(Tc_init_sim.t,Tc_init_sim.x)
plt.grid()
plt.ylabel('Cooling temperature')
plt.xlabel('time')
plt.show()

Look at the plots and make sure you understand the effect of the filter. Think about alternative, better ways to chose the input profile. Also, try to increase the value 35 that was added to the target input: how much can you increase this value without experiencing sudden ignition of the reactor?

Once the initial guess is generated, we compile the model containing the optimal control problem:

cstr = oc.compile_model("CSTR.CSTR_Opt", "CSTR.mo", target='ipopt')

We will now initialize the parameters of the model so that their values correspond to the optimization objective of transferring the system state from operating point A to operating point B. Accordingly, we set the parameters representing the initial values of the states to point A and the reference values in the cost function to point B:

cstr.set_value('Tc_ref',Tc_0_B)
cstr.set_value('c_ref',c_0_B)
cstr.set_value('T_ref',T_0_B)

cstr.set_value('cstr.c_init',c_0_A)
cstr.set_value('cstr.T_init',T_0_A)

In order to solve the optimization problem, we need to specify the mesh on which the optimization is performed. The simultaneous optimization algorithm is based on a collocation method that corresponds to a fixed step implicit Runge-Kutta scheme, where the mesh defines the length of each step. Also, the number of collocation points in each element, or step, needs to be provided. This number corresponds to the stage order of the Runge-Kutta scheme. The selection of mesh is analogous to the choice of step length in a one-step algorithm for solving differential equations. Accordingly, the mesh needs to be fine-grained enough to ensure sufficiently accurate approximation of the differential constraint. For an overview of simultaneous optimization algorithms, see [2].

Collocation-based optimization algorithms often require a good initial guess in order to achieve fast convergence. Also, if the problem is non-convex, initialization is even more critical. Initial guesses can be provided in Optimica by the 'initialGuess' attribute, see the CSTR.mo file for an example for this. Notice that initialization in the case of collocation-based optimization methods means initialization of all the control and state profiles as a function of time. In some cases, it is sufficient to use constant profiles. For this purpose, the 'initialGuess' attribute works well. In more difficult cases, however, it may be necessary to initialize the profiles using simulation data, where an initial guess for the input(s) has been used to generate the profiles for the dependent variables. This approach for initializing the optimization problem is used in this tutorial.

We are now ready to solve the actual optimization problem. This is done by invoking the method optimize:

# Initialize the mesh
n_e = 100 # Number of elements 
hs = N.ones(n_e)*1./n_e # Equidistant points
n_cp = 3; # Number of collocation points in each element

(cstr,res) = optimize(cstr,alg_args={'n_e':n_e,'hs':hs,'n_cp':n_cp,'init_traj':res})

You should see the output of Ipopt in the Python shell as the algorithm iterates to find the optimal solution. Ipopt should terminate with a message like 'Optimal solution found' or 'Solved to an acceptable level' in order for an optimum to be found. Again, the arguments to the algorithm (number of elements, number of collocation points, element length vector and initial guess object) are given in a Python dictionary. The optimization result is returned in the output argument 'res'.

We can now retrieve the trajectories of the variables that we intend to plot:

# Extract variable profiles
c_res=res.get_variable_data('cstr.c')
T_res=res.get_variable_data('cstr.T')
Tc_res=res.get_variable_data('cstr.Tc')

c_ref=res.get_variable_data('c_ref') 
T_ref=res.get_variable_data('T_ref')
Tc_ref=res.get_variable_data('Tc_ref')

Finally, we plot the result using the functions available in matplotlib:

plt.figure(1)
plt.clf()
plt.hold(True)
plt.subplot(311)
plt.plot(c_res.t,c_res.x)
plt.plot(c_ref.t,c_ref.x,'--')
plt.grid()
plt.ylabel('Concentration')

plt.subplot(312)
plt.plot(T_res.t,T_res.x)
plt.plot(T_ref.t,T_ref.x,'--')
plt.grid()
plt.ylabel('Temperature')

plt.subplot(313)
plt.plot(Tc_res.t,Tc_res.x)
plt.plot(Tc_ref.t,Tc_ref.x,'--')
plt.grid()
plt.ylabel('Cooling temperature')
plt.xlabel('time')
plt.show()

Your should now see a plot as the one below:

Take a minute to analyze the optimal profiles and to answer the following questions:

Solving optimal control problems by means of direct collocation implies that the differential equation is approximated by a discrete time counterpart. The accuracy of the solution is dependent on the method of collocation and the number of elements. In order to assess the accuracy of the discretization, we may simulate the system using a DAE solver using the optimal control profile as input. With this approach, the state profiles are computed with high accuracy and the result may then be compared with the profiles resulting from optimization. Notice that this procedure does not verify the optimality of the resulting optimal control profiles, but only the accuracy of the discretization of the dynamics.

The procedure for setting up and executing this simulation is similar to above:

# Simulate to verify the optimal solution
# Set up input trajectory
t = Tc_res.t 
u = Tc_res.x
u_traj = N.transpose(N.vstack((t,u)))
    
# Comile the Modelica model first to C code and
# then to a dynamic library
sim_model = oc.compile_model("CSTR.CSTR","CSTR.mo",target='ipopt')

sim_model.set_value('c_init',c_0_A)
sim_model.set_value('T_init',T_0_A)
sim_model.set_value('Tc',u[0])

(sim_model,res) = simulate(sim_model,compiler='optimica',
                           alg_args={'start_time':0.,'final_time':150.,
                                     'input_trajectory':u_traj})

Finally, we load the simulated data and plot it to compare with the optimized trajectories:

# Extract variable profiles
c_sim=res.get_variable_data('c')
T_sim=res.get_variable_data('T')
Tc_sim=res.get_variable_data('Tc')

# Plot the results
plt.figure(3)
plt.clf()
plt.hold(True)
plt.subplot(311)
plt.plot(c_res.t,c_res.x,'--')
plt.plot(c_sim.t,c_sim.x)
plt.legend(('optimized','simulated'))
plt.grid()
plt.ylabel('Concentration')

plt.subplot(312)
plt.plot(T_res.t,T_res.x,'--')
plt.plot(T_sim.t,T_sim.x)
plt.legend(('optimized','simulated'))
plt.grid()
plt.ylabel('Temperature')

plt.subplot(313)
plt.plot(Tc_res.t,Tc_res.x,'--')
plt.plot(Tc_sim.t,Tc_sim.x)
plt.legend(('optimized','simulated'))
plt.grid()
plt.ylabel('Cooling temperature')
plt.xlabel('time')
plt.show()

You should now see a plot similar to:

Discuss why the simulated trajectories differs from the optimized counterparts.

After completing the tutorial you may continue to modify the optimization problem and study the results.

[1] G.A. Hicks and W.H. Ray. Approximation Methods for Optimal Control Synthesis. Can. J. Chem. Eng., 40:522–529, 1971.

[2] Bieger, L., A. Cervantes, and A. Wächter (2002): "Advances in simultaneous strategies for dynamic optimization." Chemical Engineering Science, 57, pp. 575-593.