source: branches/dev-5819/Python/src/pyjmi/examples/cstr_mpc_casadi.py @ 13800

Last change on this file since 13800 was 13800, checked in by randersson, 8 weeks ago

#5819 Merged trunk into branch

File size: 9.4 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3
4# Copyright (C) 2014 Modelon AB
5#
6# This program is free software: you can redistribute it and/or modify
7# it under the terms of the GNU General Public License as published by
8# the Free Software Foundation, version 3 of the License.
9#
10# This program is distributed in the hope that it will be useful,
11# but WITHOUT ANY WARRANTY; without even the implied warranty of
12# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13# GNU General Public License for more details.
14#
15# You should have received a copy of the GNU General Public License
16# along with this program. If not, see <http://www.gnu.org/licenses/>.
17
18# Import library for path manipulations
19import os.path
20
21# Import numerical libraries
22import numpy as N
23import matplotlib.pyplot as plt
24
25# Import the needed JModelica.org Python methods
26from pymodelica import compile_fmu
27from pyfmi import load_fmu
28from pyjmi import transfer_optimization_problem, get_files_path
29from pyjmi.optimization.mpc import MPC
30from pyjmi.optimization.casadi_collocation import BlockingFactors
31import copy
32
33def run_demo(with_plots=True):
34    """
35    This example is based on the Hicks-Ray Continuously Stirred Tank Reactors
36    (CSTR) system. The system has two states, the concentration and the
37    temperature. The control input to the system is the temperature of the
38    cooling flow in the reactor jacket. The chemical reaction in the reactor is
39    exothermic, and also temperature dependent; high temperature results in
40    high reaction rate.
41
42    The problem is solved using the CasADi-based collocation algorithm through
43    the MPC-class. FMI is used for initialization and simulation purposes.
44
45    The following steps are demonstrated in this example:
46
47    1.  How to generate an initial guess for a direct collocation method by
48        means of simulation with a constant input. The trajectories resulting
49        from the simulation are used to initialize the variables in the
50        transcribed NLP, in the first sample(optimization).
51
52    2.  An optimal control problem is defined where the objective is to
53        transfer the state of the system from stationary point A to point B.
54        An MPC object for the optimization problem is created. After each
55        sample the NLP is updated with an estimate of the states in the next
56        sample. The estimate is done by simulating the model for one sample
57        period with the optimal input calculated in the optimization as input.
58        To each estimate a normally distributed noise, with the mean 0
59        and standard deviation 0.5% of the nominal value of each state, is
60        added. The MPC object uses the result from the previous optimization as
61        initial guess for the next optimization (for all but the first
62        optimization, where the simulation result from #1 is used instead).
63
64   (3.) If with_plots is True we compile the same optimization problem again
65        and define the options so that the op has the same options and
66        resolution as the op we solved through the MPC-class. By same
67        resolution we mean that both op should have the same mesh and blocking
68        factors. This allows us to compare the MPC-results to an open loop
69        optimization. Note that the MPC-results contains noise while the open
70        loop optimization does not.
71
72    """
73    ### 1. Compute initial guess trajectories by means of simulation
74    # Locate the Modelica and Optimica code
75    file_path = os.path.join(get_files_path(), "CSTR.mop")
76
77    # Compile and load the model used for simulation
78    sim_fmu = compile_fmu("CSTR.CSTR_MPC_Model", file_path, 
79                            compiler_options={"state_initial_equations":True})
80    sim_model = load_fmu(sim_fmu)
81
82    # Define stationary point A and set initial values and inputs
83    c_0_A = 956.271352
84    T_0_A = 250.051971
85    sim_model.set('_start_c', c_0_A)
86    sim_model.set('_start_T', T_0_A)
87    sim_model.set('Tc', 280)
88   
89    opts = sim_model.simulate_options()
90    opts["CVode_options"]["maxh"] = 0.0
91    opts["ncp"] = 0
92   
93    init_res = sim_model.simulate(start_time=0., final_time=150, options=opts)
94
95    ### 2. Define the optimal control problem and solve it using the MPC class
96    # Compile and load optimization problem
97    op = transfer_optimization_problem("CSTR.CSTR_MPC", file_path,
98                            compiler_options={"state_initial_equations":True})
99
100    # Define MPC options
101    sample_period = 3                           # s
102    horizon = 33                                # Samples on the horizon
103    n_e_per_sample = 1                          # Collocation elements / sample
104    n_e = n_e_per_sample*horizon                # Total collocation elements
105    finalTime = 150                             # s
106    number_samp_tot = int(finalTime/sample_period)   # Total number of samples to do
107
108    # Create blocking factors with quadratic penalty and bound on 'Tc'
109    bf_list = [n_e_per_sample]*(horizon/n_e_per_sample)
110    factors = {'Tc': bf_list}
111    du_quad_pen = {'Tc': 500}
112    du_bounds = {'Tc': 30}
113    bf = BlockingFactors(factors, du_bounds, du_quad_pen)
114
115    # Set collocation options
116    opt_opts = op.optimize_options()
117    opt_opts['n_e'] = n_e
118    opt_opts['n_cp'] = 2
119    opt_opts['init_traj'] = init_res
120    opt_opts['blocking_factors'] = bf
121
122    if with_plots:
123        # Compile and load a new instance of the op to compare the MPC results
124        # with an open loop optimization
125        op_open_loop = transfer_optimization_problem(
126            "CSTR.CSTR_MPC", file_path,
127            compiler_options={"state_initial_equations":True})
128        op_open_loop.set('_start_c', float(c_0_A))
129        op_open_loop.set('_start_T', float(T_0_A)) 
130       
131        # Copy options from MPC optimization
132        open_loop_opts = copy.deepcopy(opt_opts)
133       
134        # Change n_e and blocking_factors so op_open_loop gets the same
135        # resolution as op
136        open_loop_opts['n_e'] = number_samp_tot
137       
138        bf_list_ol = [n_e_per_sample]*(number_samp_tot/n_e_per_sample)
139        factors_ol = {'Tc': bf_list_ol}
140        bf_ol = BlockingFactors(factors_ol, du_bounds, du_quad_pen)
141        open_loop_opts['blocking_factors'] = bf_ol
142        open_loop_opts['IPOPT_options']['print_level'] = 0
143
144    constr_viol_costs = {'T': 1e6}
145
146    # Create the MPC object
147    MPC_object = MPC(op, opt_opts, sample_period, horizon, 
148                    constr_viol_costs=constr_viol_costs, noise_seed=1)
149
150    # Set initial state
151    x_k = {'_start_c': c_0_A, '_start_T': T_0_A }
152
153    # Update the state and optimize number_samp_tot times
154    for k in range(number_samp_tot):
155
156        # Update the state and compute the optimal input for next sample period
157        MPC_object.update_state(x_k)
158        u_k = MPC_object.sample()
159
160        # Reset the model and set the new initial states before simulating
161        # the next sample period with the optimal input u_k
162        sim_model.reset()
163        sim_model.set(list(x_k.keys()), list(x_k.values()))
164        sim_res = sim_model.simulate(start_time=k*sample_period, 
165                                     final_time=(k+1)*sample_period, 
166                                     input=u_k, options=opts)
167
168        # Extract state at end of sample_period from sim_res and add Gaussian
169        # noise with mean 0 and standard deviation 0.005*(state_current_value)
170        x_k = MPC_object.extract_states(sim_res, mean=0, st_dev=0.005)
171
172
173    # Extract variable profiles
174    MPC_object.print_solver_stats()
175    complete_result = MPC_object.get_complete_results()
176    c_res_comp = complete_result['c']
177    T_res_comp = complete_result['T']
178    Tc_res_comp = complete_result['Tc']
179    time_res_comp = complete_result['time']
180
181    # Verify solution for testing purposes
182    try:
183        import casadi
184    except:
185        pass
186    else:
187        Tc_norm = N.linalg.norm(Tc_res_comp) / N.sqrt(len(Tc_res_comp))
188        assert(N.abs(Tc_norm - 311.7362) < 1e-3)
189        c_norm = N.linalg.norm(c_res_comp) / N.sqrt(len(c_res_comp))
190        assert(N.abs(c_norm - 653.5369) < 1e-3)
191        T_norm = N.linalg.norm(T_res_comp) / N.sqrt(len(T_res_comp))
192        assert(N.abs(T_norm - 328.0852) < 1e-3)
193   
194    # Plot the results
195    if with_plots: 
196        ### 3. Solve the original optimal control problem without MPC
197        res = op_open_loop.optimize(options=open_loop_opts)
198        c_res = res['c']
199        T_res = res['T']
200        Tc_res = res['Tc']
201        time_res = res['time']
202       
203        # Get reference values
204        Tc_ref = op.get('Tc_ref')
205        T_ref = op.get('T_ref')
206        c_ref = op.get('c_ref')
207
208        # Plot
209        plt.close('MPC')
210        plt.figure('MPC')
211        plt.subplot(3, 1, 1)
212        plt.plot(time_res_comp, c_res_comp)
213        plt.plot(time_res, c_res )
214        plt.plot([time_res[0],time_res[-1]],[c_ref,c_ref],'--')
215        plt.legend(('MPC with noise', 'Open-loop without noise', 'Reference value'))
216        plt.grid()
217        plt.ylabel('Concentration')
218        plt.title('Simulated trajectories')
219
220        plt.subplot(3, 1, 2)
221        plt.plot(time_res_comp, T_res_comp)
222        plt.plot(time_res, T_res)
223        plt.plot([time_res[0],time_res[-1]],[T_ref,T_ref], '--')
224        plt.grid()
225        plt.ylabel('Temperature [C]')
226
227        plt.subplot(3, 1, 3)
228        plt.step(time_res_comp, Tc_res_comp)
229        plt.step(time_res, Tc_res)
230        plt.plot([time_res[0],time_res[-1]],[Tc_ref,Tc_ref], '--')
231        plt.grid()
232        plt.ylabel('Cooling temperature [C]')
233        plt.xlabel('time')
234        plt.show()
235       
236
237if __name__=="__main__":
238    run_demo()
Note: See TracBrowser for help on using the repository browser.