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 |
---|
19 | import os.path |
---|
20 | |
---|
21 | # Import numerical libraries |
---|
22 | import numpy as N |
---|
23 | import matplotlib.pyplot as plt |
---|
24 | |
---|
25 | # Import the needed JModelica.org Python methods |
---|
26 | from pymodelica import compile_fmu |
---|
27 | from pyfmi import load_fmu |
---|
28 | from pyjmi import transfer_optimization_problem, get_files_path |
---|
29 | from pyjmi.optimization.mpc import MPC |
---|
30 | from pyjmi.optimization.casadi_collocation import BlockingFactors |
---|
31 | import copy |
---|
32 | |
---|
33 | def 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 | |
---|
237 | if __name__=="__main__": |
---|
238 | run_demo() |
---|