ocp/hydroscal.cpp
Go to the documentation of this file.
1 /*
2  * This file is part of ACADO Toolkit.
3  *
4  * ACADO Toolkit -- A Toolkit for Automatic Control and Dynamic Optimization.
5  * Copyright (C) 2008-2014 by Boris Houska, Hans Joachim Ferreau,
6  * Milan Vukov, Rien Quirynen, KU Leuven.
7  * Developed within the Optimization in Engineering Center (OPTEC)
8  * under supervision of Moritz Diehl. All rights reserved.
9  *
10  * ACADO Toolkit is free software; you can redistribute it and/or
11  * modify it under the terms of the GNU Lesser General Public
12  * License as published by the Free Software Foundation; either
13  * version 3 of the License, or (at your option) any later version.
14  *
15  * ACADO Toolkit is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  * Lesser General Public License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with ACADO Toolkit; if not, write to the Free Software
22  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23  *
24  */
25 
26 
27 
54 #include <acado_optimal_control.hpp>
55 #include <acado_gnuplot.hpp>
56 #include "../integrator/hydroscal_model.hpp" // MODEL FILE
57 
58 #include <time.h>
59 
60 USING_NAMESPACE_ACADO
61 
62 void ffcn_model( double *x, double *f, void *user_data ){
63 
64  int i;
65  double *xd = new double[NXD];
66  double *xa = new double[NXA];
67  double *u = new double[ NU];
68  double *p = new double[ NP];
69 
70  for( i = 0; i < NXD; i++ ) xd[i] = x[ 1+i ];
71  for( i = 0; i < NXA; i++ ) xa[i] = x[ NXD+1+i ];
72  for( i = 0; i < NU; i++ ) u[i] = x[ NXA+NXD+1+i ];
73  for( i = 0; i < NP; i++ ) p[i] = x[ NXA+NXD+NU+1+i ];
74 
75  ffcn( &x[0], xd, xa, u, p, f );
76  gfcn( &x[0], xd, xa, u, p, &(f[NXD]) );
77 
78  delete[] xd;
79  delete[] xa;
80  delete[] u;
81  delete[] p;
82 
83 }
84 
85 
86 
87 int main( ){
88 
89  double clock1 = clock();
90 
91  double t_start = 0.0;
92  double t_end = 600.0;
93  int intervals = 5 ;
94  int i;
95 
96  TIME t;
97  DifferentialState x("", NXD, 1);
98  AlgebraicState z("", NXA, 1);
99  Control u("", NU, 1);
100  Parameter p("", NP, 1);
101  IntermediateState is(1+NXD+NXA+NU+NP);
102 
103  is(0) = t;
104  for (i=0; i < NXD; ++i) is(1+i) = x(i);
105  for (i=0; i < NXA; ++i) is(1+NXD+i) = z(i);
106  for (i=0; i < NU; ++i) is(1+NXD+NXA+i) = u(i);
107  for (i=0; i < NP; ++i) is(1+NXD+NXA+NU+i) = p(i);
108 
109  CFunction hydroscalModel( NXD+NXA, ffcn_model );
110 
111 
112  // Define a Right-Hand-Side:
113  // -------------------------
114  DifferentialEquation f;
115  f << hydroscalModel(is);
116 
117 
118  // DEFINE INITIAL VALUES:
119  // ----------------------
120  double xd[NXD] = { 2.1936116177990631E-01, // X_0
121  3.3363028623863722E-01, // X_1
122  3.7313133250625952E-01, // X_2
123  3.9896472354654333E-01, // X_3
124  4.1533719381260475E-01, // X_4
125  4.2548399372287182E-01, // X_5
126 
127  4.3168379354213621E-01, // X_6
128  4.3543569751236455E-01, // X_7
129  4.3768918647214428E-01, // X_8
130  4.3903262905928286E-01, // X_9
131  4.3982597315656735E-01, // X_10
132  4.4028774979047969E-01, // X_11
133  4.4055002518902953E-01, // X_12
134  4.4069238917008052E-01, // X_13
135  4.4076272408112094E-01, // X_14
136  4.4078980543461005E-01, // X_15
137  4.4079091412311144E-01, // X_16
138  4.4077642312834125E-01, // X_17
139  4.4075255679998443E-01, // X_18
140  4.4072304911231042E-01, // X_19
141  4.4069013958173919E-01, // X_20
142  6.7041926189645151E-01, // X_21
143  7.3517997375758948E-01, // X_22
144  7.8975978943631409E-01, // X_23
145  8.3481725159539033E-01, // X_24
146  8.7125377077380739E-01, // X_25
147  9.0027275078767721E-01, // X_26
148  9.2312464536394301E-01, // X_27
149  9.4096954980798608E-01, // X_28
150  9.5481731262797742E-01, // X_29
151  9.6551271145368878E-01, // X_30
152  9.7374401773010488E-01, // X_31
153  9.8006186072166701E-01, // X_32
154  9.8490109485675337E-01, // X_33
155  9.8860194771099286E-01, // X_34
156  9.9142879342008328E-01, // X_35
157  9.9358602331847468E-01, // X_36
158  9.9523105632238640E-01, // X_37
159  9.9648478785701988E-01, // X_38
160  9.9743986301741971E-01, // X_39
161  9.9816716097314861E-01, // X_40
162  9.9872084014280071E-01, // X_41
163  3.8633811956730968E+00, // n_1
164  3.9322260498028840E+00, // n_2
165  3.9771965626392531E+00, // n_3
166  4.0063070333869728E+00, // n_4
167  4.0246026844143410E+00, // n_5
168  4.0358888958821835E+00, // n_6
169  4.0427690398786789E+00, // n_7
170  4.0469300433477020E+00, // n_8
171  4.0494314648020326E+00, // n_9
172  4.0509267560029381E+00, // n_10
173  4.0518145583397631E+00, // n_11
174  4.0523364846379799E+00, // n_12
175  4.0526383977460299E+00, // n_13
176  4.0528081437632766E+00, // n_14
177  4.0528985491134542E+00, // n_15
178  4.0529413510270169E+00, // n_16
179  4.0529556049324462E+00, // n_17
180  4.0529527471448805E+00, // n_18
181  4.0529396392278008E+00, // n_19
182  4.0529203970496912E+00, // n_20
183  3.6071164950918582E+00, // n_21
184  3.7583754503438387E+00, // n_22
185  3.8917148481441974E+00, // n_23
186  4.0094300698741563E+00, // n_24
187  4.1102216725798293E+00, // n_25
188  4.1944038520620675E+00, // n_26
189  4.2633275166560596E+00, // n_27
190  4.3188755452109175E+00, // n_28
191  4.3630947909857642E+00, // n_29
192  4.3979622247841386E+00, // n_30
193  4.4252580012497740E+00, // n_31
194  4.4465128947193868E+00, // n_32
195  4.4630018314791968E+00, // n_33
196  4.4757626150015568E+00, // n_34
197  4.4856260094946823E+00, // n_35
198  4.4932488551808500E+00, // n_36
199  4.4991456959629330E+00, // n_37
200  4.5037168116896273E+00, // n_38
201  4.5072719605639726E+00, // n_39
202  4.5100498969782414E+00 // n_40
203  };
204 
205  double ud[ NU] = { 4.1833910982822058E+00, // L_vol
206  2.4899344742988991E+00 // Q
207  };
208 
209  double pd[ NP] = { 1.5458567140000001E-01, // n^{ref}_{tray} = 0.155 l
210  1.7499999999999999E-01, // (not in use?)
211  3.4717208398678062E-01, // \alpha_{rect} = 35 %
212  6.1895708603484367E-01, // \alpha_{strip} = 62 %
213  1.6593025789999999E-01, // W_{tray} = 0.166 l^{-.5} s^{.1}
214  5.0695122527590109E-01, // Q_{loss} = 0.51 kW
215  8.5000000000000000E+00, // n^v_0 = 8.5 l
216  1.7000000000000001E-01, // n^v_{N+1} = 0.17 l
217  9.3885430857029321E+04, // P_{top} = 939 h Pa
218  2.5000000000000000E+02, // \Delta P_{strip}= 2.5 h Pa and \Delta P_{rect} = 1.9 h Pa
219  1.4026000000000000E+01, // F_{vol} = 14.0 l h^{-1}
220  3.2000000000000001E-01, // X_F = 0.32
221  7.1054000000000002E+01, // T_F = 71 oC
222  4.7163089489100003E+01, // T_C = 47.2 oC
223  4.1833910753991770E+00, // (not in use?)
224  2.4899344810136301E+00, // (not in use?)
225  1.8760537088149468E+02 // (not in use?)
226  };
227 
228  DVector x0(NXD, xd);
229  DVector p0(NP, pd);
230 
231 
232  // DEFINE AN OPTIMAL CONTROL PROBLEM:
233  // ----------------------------------
234  OCP ocp( t_start, t_end, intervals );
235 
236  // LSQ Term on temperature deviations and controls
237  Function h;
238 
239 
240 // for( i = 0; i < NXD; i++ )
241 // h << 0.001*x(i);
242 
243 
244  h << 0.1 * ( z(94) - 88.0 ); // Temperature tray 14
245  h << 0.1 * ( z(108) - 70.0 ); // Temperature tray 28
246  h << 0.01 * ( u(0) - ud[0] ); // L_vol
247  h << 0.01 * ( u(1) - ud[1] ); // Q
248  ocp.minimizeLSQ( h );
249 
250  // W.r.t. differential equation
251  ocp.subjectTo( f );
252 
253  // Fix states
254  ocp.subjectTo( AT_START, x == x0 );
255 
256  // Fix parameters
257  ocp.subjectTo( p == p0 );
258 
259  // Path constraint on controls
260  ocp.subjectTo( ud[0] - 2.0 <= u(0) <= ud[0] + 2.0 );
261  ocp.subjectTo( ud[1] - 2.0 <= u(1) <= ud[1] + 2.0 );
262 
263 
264 
265  // DEFINE AN OPTIMIZATION ALGORITHM AND SOLVE THE OCP:
266  // ---------------------------------------------------
267  OptimizationAlgorithm algorithm(ocp);
268 
269  algorithm.initializeAlgebraicStates("hydroscal_algebraic_states.txt");
270 
271  algorithm.set( INTEGRATOR_TYPE, INT_BDF );
272  algorithm.set( MAX_NUM_ITERATIONS, 5 );
273  algorithm.set( KKT_TOLERANCE, 1e-3 );
274  algorithm.set( INTEGRATOR_TOLERANCE, 1e-4 );
275  algorithm.set( ABSOLUTE_TOLERANCE , 1e-6 );
276  algorithm.set( PRINT_SCP_METHOD_PROFILE, YES );
277  algorithm.set( LINEAR_ALGEBRA_SOLVER, SPARSE_LU );
278  algorithm.set( DISCRETIZATION_TYPE, MULTIPLE_SHOOTING );
279 
280  //algorithm.set( LEVENBERG_MARQUARDT, 1e-3 );
281 
282  algorithm.set( DYNAMIC_SENSITIVITY, FORWARD_SENSITIVITY_LIFTED );
283  //algorithm.set( DYNAMIC_SENSITIVITY, FORWARD_SENSITIVITY );
284  //algorithm.set( CONSTRAINT_SENSITIVITY, FORWARD_SENSITIVITY );
285  //algorithm.set( ALGEBRAIC_RELAXATION,ART_EXPONENTIAL ); //results in an extra step but steps are quicker
286 
287 
288  algorithm.solve();
289 
290  double clock2 = clock();
291  printf("total computation time = %.16e \n", (clock2-clock1)/CLOCKS_PER_SEC );
292 
293 
294  // PLOT RESULTS:
295  // ---------------------------------------------------
296  VariablesGrid out_states;
297  algorithm.getDifferentialStates( out_states );
298  out_states.print( "OUT_states.m","STATES",PS_MATLAB );
299 
300  VariablesGrid out_controls;
301  algorithm.getControls( out_controls );
302  out_controls.print( "OUT_controls.m","CONTROLS",PS_MATLAB );
303 
304  VariablesGrid out_algstates;
305  algorithm.getAlgebraicStates( out_algstates );
306  out_algstates.print( "OUT_algstates.m","ALGSTATES",PS_MATLAB );
307 
308  GnuplotWindow window;
309  window.addSubplot( out_algstates(94), "Temperature tray 14" );
310  window.addSubplot( out_algstates(108), "Temperature tray 28" );
311  window.addSubplot( out_controls(0), "L_vol" );
312  window.addSubplot( out_controls(1), "Q" );
313  window.plot( );
314 
315  return 0;
316 }
317 
318 
319 
NXD
#define NXD
Definition: hydroscal_model.hpp:72
MatrixVariablesGrid::print
returnValue print(std::ostream &stream=std::cout, const char *const name=DEFAULT_LABEL, const char *const startString=DEFAULT_START_STRING, const char *const endString=DEFAULT_END_STRING, uint width=DEFAULT_WIDTH, uint precision=DEFAULT_PRECISION, const char *const colSeparator=DEFAULT_COL_SEPARATOR, const char *const rowSeparator=DEFAULT_ROW_SEPARATOR) const
Definition: matrix_variables_grid.cpp:695
ffcn_model
USING_NAMESPACE_ACADO void ffcn_model(double *x, double *f, void *user_data)
Definition: ocp/hydroscal.cpp:62
Function
Allows to setup and evaluate a general function based on SymbolicExpressions.
Definition: function_.hpp:59
PlotWindow::plot
virtual returnValue plot(PlotFrequency _frequency=PLOT_IN_ANY_CASE)
Definition: plot_window.cpp:166
acado_optimal_control.hpp
FORWARD_SENSITIVITY_LIFTED
Definition: acado_types.hpp:587
OptimizationAlgorithm
User-interface to formulate and solve optimal control problems and static NLPs.
Definition: optimization_algorithm.hpp:54
USING_NAMESPACE_ACADO
#define USING_NAMESPACE_ACADO
Definition: acado_namespace_macros.hpp:59
VariablesGrid
Provides a time grid consisting of vector-valued optimization variables at each grid point...
Definition: variables_grid.hpp:56
ABSOLUTE_TOLERANCE
Definition: acado_types.hpp:358
PS_MATLAB
Definition: acado_types.hpp:578
PRINT_SCP_METHOD_PROFILE
Definition: acado_types.hpp:418
OCP::subjectTo
returnValue subjectTo(const DifferentialEquation &differentialEquation_)
Definition: ocp.cpp:153
PlotWindow::addSubplot
returnValue addSubplot(PlotWindowSubplot &_subplot)
Definition: plot_window.cpp:276
DISCRETIZATION_TYPE
Definition: acado_types.hpp:383
OptimizationAlgorithmBase::getAlgebraicStates
returnValue getAlgebraicStates(VariablesGrid &xa_) const
Definition: optimization_algorithm_base.cpp:232
Options::set
returnValue set(OptionsName name, int value)
Definition: options.cpp:126
OCP::minimizeLSQ
returnValue minimizeLSQ(const DMatrix &S, const Function &h, const DVector &r)
Definition: ocp.cpp:244
DYNAMIC_SENSITIVITY
Definition: acado_types.hpp:380
YES
#define YES
Definition: acado_types.hpp:51
INTEGRATOR_TOLERANCE
Definition: acado_types.hpp:355
NXA
#define NXA
Definition: crane_kul_mhe_test.cpp:42
gfcn
static void gfcn(double *t, double *xd, double *xa, double *u, double *p, double *rhs)
Definition: hydroscal_model.hpp:542
OCP
Data class for defining optimal control problems.
Definition: ocp.hpp:89
SPARSE_LU
Definition: acado_types.hpp:152
ffcn
static void ffcn(double *t, double *xd, double *xa, double *u, double *p, double *rhs)
Definition: hydroscal_model.hpp:502
OptimizationAlgorithmBase::getControls
returnValue getControls(VariablesGrid &p_) const
Definition: optimization_algorithm_base.cpp:262
OptimizationAlgorithmBase::initializeAlgebraicStates
returnValue initializeAlgebraicStates(const char *fileName, BooleanType autoinit=BT_FALSE)
Definition: optimization_algorithm_base.cpp:118
KKT_TOLERANCE
Definition: acado_types.hpp:371
t_end
const double t_end
Definition: rocket_with_templates_ocp.hpp:53
MULTIPLE_SHOOTING
Definition: acado_types.hpp:203
t_start
const double t_start
Definition: rocket_with_templates_ocp.hpp:52
NP
#define NP
Definition: hydroscal_model.hpp:62
CFunction
(no description yet)
Definition: c_function.hpp:54
NU
#define NU
Definition: crane_kul_mhe_test.cpp:43
LINEAR_ALGEBRA_SOLVER
Definition: acado_types.hpp:365
MAX_NUM_ITERATIONS
Definition: acado_types.hpp:370
INT_BDF
Definition: acado_types.hpp:262
GnuplotWindow
Provides an interface to Gnuplot for plotting algorithmic outputs.
Definition: gnuplot_window.hpp:58
AT_START
Definition: acado_types.hpp:732
main
int main()
Definition: ocp/hydroscal.cpp:87
INTEGRATOR_TYPE
Definition: acado_types.hpp:397
OptimizationAlgorithm::solve
virtual returnValue solve()
Definition: optimization_algorithm.cpp:100
acado_gnuplot.hpp
OptimizationAlgorithmBase::getDifferentialStates
returnValue getDifferentialStates(VariablesGrid &xd_) const
Definition: optimization_algorithm_base.cpp:225
DifferentialEquation
Allows to setup and evaluate differential equations (ODEs and DAEs) based on SymbolicExpressions.
Definition: differential_equation.hpp:55

acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:34:41