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

Last change on this file since 13461 was 13461, checked in by randersson, 3 months ago

#5819 Performed 2to3 on the entire src/Python directory in order to make the code compatible. I've not reviewed all changes while doing this commit but will do while starting to run the tests. I'm not sure everything is converted correctly, for example the change in test_casadi_collocation row 1610.

File size: 9.3 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    init_res = sim_model.simulate(start_time=0., final_time=150)
89
90    ### 2. Define the optimal control problem and solve it using the MPC class
91    # Compile and load optimization problem
92    op = transfer_optimization_problem("CSTR.CSTR_MPC", file_path,
93                            compiler_options={"state_initial_equations":True})
94
95    # Define MPC options
96    sample_period = 3                           # s
97    horizon = 33                                # Samples on the horizon
98    n_e_per_sample = 1                          # Collocation elements / sample
99    n_e = n_e_per_sample*horizon                # Total collocation elements
100    finalTime = 150                             # s
101    number_samp_tot = int(finalTime/sample_period)   # Total number of samples to do
102
103    # Create blocking factors with quadratic penalty and bound on 'Tc'
104    bf_list = [n_e_per_sample]*(horizon/n_e_per_sample)
105    factors = {'Tc': bf_list}
106    du_quad_pen = {'Tc': 500}
107    du_bounds = {'Tc': 30}
108    bf = BlockingFactors(factors, du_bounds, du_quad_pen)
109
110    # Set collocation options
111    opt_opts = op.optimize_options()
112    opt_opts['n_e'] = n_e
113    opt_opts['n_cp'] = 2
114    opt_opts['init_traj'] = init_res
115    opt_opts['blocking_factors'] = bf
116
117    if with_plots:
118        # Compile and load a new instance of the op to compare the MPC results
119        # with an open loop optimization
120        op_open_loop = transfer_optimization_problem(
121            "CSTR.CSTR_MPC", file_path,
122            compiler_options={"state_initial_equations":True})
123        op_open_loop.set('_start_c', float(c_0_A))
124        op_open_loop.set('_start_T', float(T_0_A)) 
125       
126        # Copy options from MPC optimization
127        open_loop_opts = copy.deepcopy(opt_opts)
128       
129        # Change n_e and blocking_factors so op_open_loop gets the same
130        # resolution as op
131        open_loop_opts['n_e'] = number_samp_tot
132       
133        bf_list_ol = [n_e_per_sample]*(number_samp_tot/n_e_per_sample)
134        factors_ol = {'Tc': bf_list_ol}
135        bf_ol = BlockingFactors(factors_ol, du_bounds, du_quad_pen)
136        open_loop_opts['blocking_factors'] = bf_ol
137        open_loop_opts['IPOPT_options']['print_level'] = 0
138
139    constr_viol_costs = {'T': 1e6}
140
141    # Create the MPC object
142    MPC_object = MPC(op, opt_opts, sample_period, horizon, 
143                    constr_viol_costs=constr_viol_costs, noise_seed=1)
144
145    # Set initial state
146    x_k = {'_start_c': c_0_A, '_start_T': T_0_A }
147
148    # Update the state and optimize number_samp_tot times
149    for k in range(number_samp_tot):
150
151        # Update the state and compute the optimal input for next sample period
152        MPC_object.update_state(x_k)
153        u_k = MPC_object.sample()
154
155        # Reset the model and set the new initial states before simulating
156        # the next sample period with the optimal input u_k
157        sim_model.reset()
158        sim_model.set(list(x_k.keys()), list(x_k.values()))
159        sim_res = sim_model.simulate(start_time=k*sample_period, 
160                                     final_time=(k+1)*sample_period, 
161                                     input=u_k)
162
163        # Extract state at end of sample_period from sim_res and add Gaussian
164        # noise with mean 0 and standard deviation 0.005*(state_current_value)
165        x_k = MPC_object.extract_states(sim_res, mean=0, st_dev=0.005)
166
167
168    # Extract variable profiles
169    MPC_object.print_solver_stats()
170    complete_result = MPC_object.get_complete_results()
171    c_res_comp = complete_result['c']
172    T_res_comp = complete_result['T']
173    Tc_res_comp = complete_result['Tc']
174    time_res_comp = complete_result['time']
175
176    # Verify solution for testing purposes
177    try:
178        import casadi
179    except:
180        pass
181    else:
182        Tc_norm = N.linalg.norm(Tc_res_comp) / N.sqrt(len(Tc_res_comp))
183        assert(N.abs(Tc_norm - 311.7362) < 1e-3)
184        c_norm = N.linalg.norm(c_res_comp) / N.sqrt(len(c_res_comp))
185        assert(N.abs(c_norm - 653.5369) < 1e-3)
186        T_norm = N.linalg.norm(T_res_comp) / N.sqrt(len(T_res_comp))
187        assert(N.abs(T_norm - 328.0852) < 1e-3)
188   
189    # Plot the results
190    if with_plots: 
191        ### 3. Solve the original optimal control problem without MPC
192        res = op_open_loop.optimize(options=open_loop_opts)
193        c_res = res['c']
194        T_res = res['T']
195        Tc_res = res['Tc']
196        time_res = res['time']
197       
198        # Get reference values
199        Tc_ref = op.get('Tc_ref')
200        T_ref = op.get('T_ref')
201        c_ref = op.get('c_ref')
202
203        # Plot
204        plt.close('MPC')
205        plt.figure('MPC')
206        plt.subplot(3, 1, 1)
207        plt.plot(time_res_comp, c_res_comp)
208        plt.plot(time_res, c_res )
209        plt.plot([time_res[0],time_res[-1]],[c_ref,c_ref],'--')
210        plt.legend(('MPC with noise', 'Open-loop without noise', 'Reference value'))
211        plt.grid()
212        plt.ylabel('Concentration')
213        plt.title('Simulated trajectories')
214
215        plt.subplot(3, 1, 2)
216        plt.plot(time_res_comp, T_res_comp)
217        plt.plot(time_res, T_res)
218        plt.plot([time_res[0],time_res[-1]],[T_ref,T_ref], '--')
219        plt.grid()
220        plt.ylabel('Temperature [C]')
221
222        plt.subplot(3, 1, 3)
223        plt.step(time_res_comp, Tc_res_comp)
224        plt.step(time_res, Tc_res)
225        plt.plot([time_res[0],time_res[-1]],[Tc_ref,Tc_ref], '--')
226        plt.grid()
227        plt.ylabel('Cooling temperature [C]')
228        plt.xlabel('time')
229        plt.show()
230       
231
232if __name__=="__main__":
233    run_demo()
Note: See TracBrowser for help on using the repository browser.