plug_flow_reactor_ws.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 
46 #include <acado_gnuplot.hpp>
47 
48 
49 int main( ){
50 
52 
53  // INTRODUCE FIXED PARAMETERS:
54  // ---------------------------
55  #define v 0.1
56  #define L 1.0
57  #define Beta 0.2
58  #define Delta 0.25
59  #define E 11250.0
60  #define k0 1E+06
61  #define R 1.986
62  #define K1 250000.0
63  #define Cin 0.02
64  #define Tin 340.0
65 
66 
67  // INTRODUCE THE VARIABLES:
68  // -------------------------
69  DifferentialState x1,x2;
70  Control u ;
71  DifferentialEquation f( 0.0, L );
72 
73 
74  // DEFINE A DIFFERENTIAL EQUATION:
75  // -------------------------------
76  double Alpha, Gamma;
77  Alpha = k0*exp(-E/(R*Tin));
78  Gamma = E/(R*Tin);
79 
80  f << dot(x1) == Alpha /v * (1.0-x1) * exp((Gamma*x2)/(1.0+x2));
81  f << dot(x2) == (Alpha*Delta)/v * (1.0-x1) * exp((Gamma*x2)/(1.0+x2)) + Beta/v * (u-x2);
82 
83 
84  // DEFINE AN OPTIMAL CONTROL PROBLEM:
85  // ----------------------------------
86  OCP ocp( 0.0, L, 50 );
87  ocp.minimizeMayerTerm( 0, Cin*(1.0-x1) ); // Solve conversion optimal problem
88  ocp.minimizeMayerTerm( 1, (pow((Tin*x2),2.0)/K1) + 0.005*Cin*(1.0-x1) ); // Solve energy optimal problem (perturbed by small conversion cost;
89  // otherwise the problem is ill-defined.)
90 
91  ocp.subjectTo( f );
92 
93  ocp.subjectTo( AT_START, x1 == 0.0 );
94  ocp.subjectTo( AT_START, x2 == 0.0 );
95 
96  ocp.subjectTo( 0.0 <= x1 <= 1.0 );
97  ocp.subjectTo( (280.0-Tin)/Tin <= x2 <= (400.0-Tin)/Tin );
98  ocp.subjectTo( (280.0-Tin)/Tin <= u <= (400.0-Tin)/Tin );
99 
100 
101  // DEFINE A MULTI-OBJECTIVE ALGORITHM AND SOLVE THE OCP:
102  // -----------------------------------------------------
103  MultiObjectiveAlgorithm algorithm(ocp);
104 
105  algorithm.set( INTEGRATOR_TYPE, INT_BDF );
106  algorithm.set( KKT_TOLERANCE, 1e-8 );
107 
109  algorithm.set( PARETO_FRONT_DISCRETIZATION, 11 );
110 
111  // Generate Pareto set
112  algorithm.solve();
113 
114  algorithm.getWeights("plug_flow_reactor_ws_weights.txt");
115  algorithm.getAllDifferentialStates("plug_flow_reactor_ws_states.txt");
116  algorithm.getAllControls("plug_flow_reactor_ws_controls.txt");
117 
118 
119  // VISUALIZE THE RESULTS IN A GNUPLOT WINDOW:
120  // ------------------------------------------
121  VariablesGrid paretoFront;
122  algorithm.getParetoFront( paretoFront );
123 
124  GnuplotWindow window1;
125  window1.addSubplot( paretoFront, "Pareto Front (conversion versus energy)", "OUTLET CONCENTRATION", "ENERGY",PM_POINTS );
126  window1.plot( );
127 
128 
129  // PRINT INFORMATION ABOUT THE ALGORITHM:
130  // --------------------------------------
131  algorithm.printInfo();
132 
133 
134  // SAVE INFORMATION:
135  // -----------------
136  paretoFront.print();
137 
138  return 0;
139 }
140 
141 
142 
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
#define Beta
DMatrix getWeights() const
virtual returnValue plot(PlotFrequency _frequency=PLOT_IN_ANY_CASE)
#define USING_NAMESPACE_ACADO
Provides a time grid consisting of vector-valued optimization variables at each grid point...
#define Cin
#define Delta
returnValue printInfo()
IntermediateState pow(const Expression &arg1, const Expression &arg2)
#define k0
returnValue addSubplot(PlotWindowSubplot &_subplot)
#define E
int main()
returnValue set(OptionsName name, int value)
Definition: options.cpp:126
returnValue getAllControls(const char *fileName) const
returnValue getParetoFront(VariablesGrid &paretoFront) const
#define L
Data class for defining optimal control problems.
Definition: ocp.hpp:89
Expression dot(const Expression &arg)
#define Tin
#define R
#define v
IntermediateState exp(const Expression &arg)
#define K1
User-interface to formulate and solve optimal control problems with multiple objectives.
returnValue getAllDifferentialStates(const char *fileName) const
Provides an interface to Gnuplot for plotting algorithmic outputs.
Allows to setup and evaluate differential equations (ODEs and DAEs) based on SymbolicExpressions.


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