multi_objective_algorithm.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 
35 #include <acado/ocp/ocp.hpp>
36 
38 
39 
40 //
41 // PUBLIC MEMBER FUNCTIONS:
42 //
43 
44 
47 
48  m = 0;
49  N = 0;
50  count = 0;
52  totalCPUtime = 0;
53 
54  xResults = 0;
55  xaResults = 0;
56  pResults = 0;
57  uResults = 0;
58  wResults = 0;
59 
60  setupOptions( );
61 }
62 
63 
65  :OptimizationAlgorithm( ocp_ ){
66 
68  N = 0;
69  count = 0;
71  totalCPUtime = 0;
72 
73  xResults = 0;
74  xaResults = 0;
75  pResults = 0;
76  uResults = 0;
77  wResults = 0;
78 
79  setupOptions( );
80 }
81 
82 
84  :OptimizationAlgorithm( arg )
85 {
86 
87  m = arg.m;
88  N = arg.N;
89 
90  vertices = arg.vertices;
91  result = arg.result ;
92 
93  count = arg.count;
94 
97 
98  xResults = 0;
99  xaResults = 0;
100  pResults = 0;
101  uResults = 0;
102  wResults = 0;
103 }
104 
105 
107 
108  if( xResults != 0 ) delete[] xResults ;
109  if( xaResults != 0 ) delete[] xaResults;
110  if( pResults != 0 ) delete[] pResults ;
111  if( uResults != 0 ) delete[] uResults ;
112  if( wResults != 0 ) delete[] wResults ;
113 }
114 
115 
117 
118  if( this != &arg ){
119 
121 
122  m = arg.m;
123  N = arg.N;
124 
125  vertices = arg.vertices;
126  result = arg.result ;
127  count = arg.count;
128 
130  totalCPUtime = arg.totalCPUtime ;
131 
132  xResults = 0;
133  xaResults = 0;
134  pResults = 0;
135  uResults = 0;
136  wResults = 0;
137  }
138  return *this;
139 }
140 
141 
143 
144 
145  int run1 ;
146  returnValue returnvalue;
147 
148  ASSERT( ocp != 0 );
149  ASSERT( number_ < m );
150  if( N == 0 ) get( PARETO_FRONT_DISCRETIZATION, N );
151 
152  int N_pow_m_minus_1 = 1;
153  int i;
154  for(i=0; i<m-1; i++)
155  N_pow_m_minus_1 *= N;
156 
157  Expression **arg = 0;
158  arg = new Expression*[m];
159  for( run1 = 0; run1 < m; run1++ )
160  ocp->getObjective( run1, &arg[run1] );
161 
162  Grid tmp_grid;
163  ocp->getGrid( tmp_grid );
164  Objective tmp(tmp_grid);
165  tmp.addMayerTerm(*arg[number_]);
166  ocp->setObjective( tmp );
167 
169 
170  printf("\n\n Optimization of individual objective %d out of %d \n\n",number_ +1, m );
172  returnvalue = OptimizationAlgorithm::solve();
174 
175  int index=0;
176  DMatrix Weights = getWeights();
177 
178  for( run1 = 0; run1 < (int) Weights.getNumCols(); run1++ ){
179  if( fabs( Weights( number_, run1 ) - 1.0 ) < 1000.0*EPS )
180  index = run1;
181  }
182 
183  if( xResults == 0 ) xResults = new VariablesGrid[Weights.getNumCols()];
184  if( xaResults == 0 ) xaResults = new VariablesGrid[Weights.getNumCols()];
185  if( pResults == 0 ) pResults = new VariablesGrid[Weights.getNumCols()];
186  if( uResults == 0 ) uResults = new VariablesGrid[Weights.getNumCols()];
187  if( wResults == 0 ) wResults = new VariablesGrid[Weights.getNumCols()];
188 
190  getAlgebraicStates ( xaResults[index] );
191  getParameters ( pResults[index] );
192  getControls ( uResults[index] );
193  getDisturbances ( wResults[index] );
194 
195 
196  if( returnvalue != SUCCESSFUL_RETURN )
197  return ACADOERROR(returnvalue);
198 
199  if( nlpSolver != 0 )
201 
202  int hotstart;
204 
205  int *indices = new int[m];
206  for( run1 = 0; run1 < m; run1++ )
207  indices[run1] = number_;
208 
209  VariablesGrid xd_tmp, xa_tmp, p_tmp, u_tmp, w_tmp;
210 
211  if( hotstart == BT_TRUE ){
214  getParameters ( *userInit.p );
215  getControls ( *userInit.u );
217  xd_tmp = *userInit.x;
218  xa_tmp = *userInit.xa;
219  p_tmp = *userInit.p;
220  u_tmp = *userInit.u;
221  w_tmp = *userInit.w;
222  }
223  else{
224  getDifferentialStates( xd_tmp );
225  getAlgebraicStates ( xa_tmp );
226  getParameters ( p_tmp );
227  getControls ( u_tmp );
228  getDisturbances ( w_tmp );
229  }
230 
231  VariablesGrid *_xd = 0;
232  VariablesGrid *_xa = 0;
233  VariablesGrid *_p = 0;
234  VariablesGrid *_u = 0;
235  VariablesGrid *_w = 0;
236 
237  if( xd_tmp.isEmpty() == BT_FALSE ) _xd = new VariablesGrid(xd_tmp);
238  if( xa_tmp.isEmpty() == BT_FALSE ) _xa = new VariablesGrid(xa_tmp);
239  if( p_tmp.isEmpty() == BT_FALSE ) _p = new VariablesGrid(p_tmp );
240  if( u_tmp.isEmpty() == BT_FALSE ) _u = new VariablesGrid(u_tmp );
241  if( w_tmp.isEmpty() == BT_FALSE ) _w = new VariablesGrid(w_tmp );
242 
243  if( vertices.getDim() == 0 ){
244  vertices.init(m,m);
246  }
247 
248  Objective **obj = new Objective*[m];
249  for( run1 = 0; run1 < m; run1++ ){
250  obj[run1] = new Objective( tmp_grid );
251  obj[run1]->addMayerTerm(*arg[run1]);
252  OCPiterate xx( _xd, _xa, _p, _u, _w );
253  obj[run1]->evaluate( xx );
254  obj[run1]->getObjectiveValue( vertices(number_,run1) );
255  }
256 
257  if( _xd != 0 ) delete _xd;
258  if( _xa != 0 ) delete _xa;
259  if( _p != 0 ) delete _p ;
260  if( _u != 0 ) delete _u ;
261  if( _w != 0 ) delete _w ;
262 
263 
264  delete[] indices;
265 
266  for( run1 = 0; run1 < m; run1++ )
267  delete arg[run1];
268  delete[] arg;
269 
270  for( run1 = 0; run1 < m; run1++ )
271  delete obj[run1];
272  delete[] obj;
273 
274  return SUCCESSFUL_RETURN;
275 }
276 
277 
279 
280  int run1,run2;
281  returnValue returnvalue;
282 
283  ASSERT( ocp != 0 );
284  ASSERT( m >= 2 );
285  if( N == 0 ) get( PARETO_FRONT_DISCRETIZATION, N );
286 
287  int paretoGeneration;
288  get( PARETO_FRONT_GENERATION, paretoGeneration );
289 
290  int hotstart;
292 
293  Expression **arg = 0;
294  arg = new Expression*[m];
295 
296  for( run1 = 0; run1 < m; run1++ )
297  ocp->getObjective( run1, &arg[run1] );
298 
299  Constraint tmp_con;
300 
301  double *idx = new double[m];
302 
303  WeightGeneration generator;
304  DMatrix Weights;
305  DVector formers;
306  DVector lb(m);
307  DVector ub(m);
308  lb.setZero();
309  ub.setAll(1.0);
310 
311  generator.getWeights( m, N, lb, ub, Weights, formers );
312 
313  result.init( Weights.getNumCols(), m );
314  count = 0;
315 
316  if( xResults == 0 ) xResults = new VariablesGrid[Weights.getNumCols()];
317  if( xaResults == 0 ) xaResults = new VariablesGrid[Weights.getNumCols()];
318  if( pResults == 0 ) pResults = new VariablesGrid[Weights.getNumCols()];
319  if( uResults == 0 ) uResults = new VariablesGrid[Weights.getNumCols()];
320  if( wResults == 0 ) wResults = new VariablesGrid[Weights.getNumCols()];
321 
324 
325  run1 = 0;
326  while( run1 < (int) Weights.getNumCols() ){
327 
328 
329  // PRINT THE ITERATION NUMBER:
330  // ---------------------------
331  printf("\n\n Multi-objective point: %d out of %d \n\n",run1+1, (int) Weights.getNumCols() );
332 
333 
334  ocp->getConstraint( tmp_con );
335 
336  for( run2 = 0; run2 < (int) Weights.getNumRows(); run2++ )
337  idx[run2] = Weights( run2, run1 );
338 
339 
340  // THIS PART OF THE CODE WILL NOT RUN YET FOR GENERAL WEIGHTS
341 
342  int vertex = -1;
343  for( run2 = 0; run2 < m; run2++ ){
344  if( fabs( idx[run2]-1.0 ) < 100.0*EPS )
345  vertex = run2;
346  }
347  // ----------------------------------------------------------
348 
349 
350  if( vertex == -1 || paretoGeneration == PFG_WEIGHTED_SUM ){
351 
352  formulateOCP( idx, ocp, arg );
354  returnvalue = OptimizationAlgorithm::solve();
355 
356  if( nlpSolver != 0 )
358 
359  ocp->setConstraint( tmp_con );
360  set( PRINT_COPYRIGHT, BT_FALSE );
361 
362  if( returnvalue != SUCCESSFUL_RETURN ){
363  ACADOERROR(returnvalue);
364  }
365  else{
366 
368  getAlgebraicStates ( xaResults[run1] );
369  getParameters ( pResults[run1] );
370  getControls ( uResults[run1] );
371  getDisturbances ( wResults[run1] );
372 
373  if( hotstart == BT_TRUE ){
376  getParameters ( *userInit.p );
377  getControls ( *userInit.u );
380  }
381  else{
382  VariablesGrid xd_tmp, xa_tmp, p_tmp, u_tmp, w_tmp;
383  getDifferentialStates( xd_tmp );
384  getAlgebraicStates ( xa_tmp );
385  getParameters ( p_tmp );
386  getControls ( u_tmp );
387  getDisturbances ( w_tmp );
388  evaluateObjectives( xd_tmp, xa_tmp, p_tmp, u_tmp, w_tmp, arg );
389  }
390  }
391  }
392  else{
393  printf(" Result from single objective optimization is adopted. \n\n" );
394  for( run2 = 0; run2 < m; run2++ ){
395  result(count,run2) = vertices(vertex,run2);
396  }
397  count++;
398  }
399  run1++;
400  }
402 
403  for( run1 = 0; run1 < m; run1++ )
404  delete arg[run1];
405  delete[] arg;
406 
407  delete[] idx;
408 
409  return SUCCESSFUL_RETURN;
410 }
411 
412 
413 
414 //
415 // PROTECTED MEMBER FUNCTIONS:
416 //
417 
419 {
423 
424  // add optimization algorithm options
425  //OptimizationAlgorithm::setupOptions( );
426 
427  return SUCCESSFUL_RETURN;
428 }
429 
430 
432 {
433  return OptimizationAlgorithm::initializeNlpSolver( _userInit );
434 }
435 
436 
438  )
439 {
440  return SUCCESSFUL_RETURN;
441 }
442 
443 
444 
445 
447 
448  int run1, run2;
449  int paretoGeneration;
450  double factor;
451 
452  get( PARETO_FRONT_GENERATION, paretoGeneration );
453  if( paretoGeneration == PFG_UNKNOWN ) paretoGeneration = PFG_WEIGHTED_SUM;
454 
455  Grid tmp_grid;
456  ocp_->getGrid( tmp_grid );
457  Objective tmp(tmp_grid);
458 
459 
460  // WEIGTHED SUM:
461  // -----------------------------------------------------------------------------
462 
463  if( paretoGeneration == PFG_WEIGHTED_SUM ){
464 
465  Expression sum(1);
466 
467  for( run1 = 0; run1 < m; run1++ ){ // loop over the number of objectives
468  factor = idx[run1]; // determines the weight factor
469  sum = sum + arg[run1][0]*factor;
470  }
471  tmp.addMayerTerm(sum); // add the new objective as a Meyer Term
472  ocp_->setObjective(tmp); // replace (overwrite) the objective in the ocp
473  return SUCCESSFUL_RETURN;
474  }
475 
476 
477  // NORMALIZED NORMAL CONSTRAINT:
478  // -----------------------------------------------------------------------------
479 
480  if( paretoGeneration == PFG_NORMALIZED_NORMAL_CONSTRAINT ){
481  // Normalization based on Messac et al 2004
482 
483  //tmp.addMayerTerm( *arg[m-1] ); // Select last (i.e., m-th) objective function
484  //ocp_->setObjective( tmp ); // replace (overwrite) the objective in the ocp
485 
488 
489  DVector W(m);
490  for( run1 = 0; run1 < m; run1++ )
491  W(run1) = idx[run1];
492 
493  DVector PW = P*W;
494 
497 
498  Expression *Fnorm;
499  Fnorm = new Expression[m];
500 
501  for( run2 = 0; run2 < m; run2++ ){
502 
503  Fnorm[run2] = ( *arg[run2] - U(run2) ) / L(run2);
504  }
505 
506  tmp.addMayerTerm( Fnorm[m-1] );
507  ocp_->setObjective( tmp );
508 
509  for( run1 = 0; run1 < m-1; run1++ ){
510  Expression sum(1);
511  for( run2 = 0; run2 < m; run2++ ){
512 
513  sum = sum + NK(run2,run1)*( Fnorm[run2] - PW(run2) );
514  }
515  ocp_->subjectTo( AT_END, sum <= 0.0 );
516  }
517  delete[] Fnorm;
518 
519  return SUCCESSFUL_RETURN;
520  }
521 
522 
523  // ENHANCED NORMALIZED NORMAL CONSTRAINT:
524  // -----------------------------------------------------------------------------
525 
526  if( paretoGeneration == PFG_ENHANCED_NORMALIZED_NORMAL_CONSTRAINT ){
527  // Normalization based on Sanchis et al 2008
528 
529  DMatrix P = getPayOffMatrix();
530  DMatrix PHI_N(m,m);
531 
532  int run3, run4;
533  for( run3 = 0; run3 < m; run3++ ){
534  for( run4 = 0; run4 < m; run4++ ){
535  PHI_N(run3,run4) = 1.0;
536  }
537  }
538  for( run3 = 0; run3 < m; run3++ )
539  PHI_N(run3,run3) = 0.0;
540 
541  DMatrix T;
542  T = PHI_N * P.inverse();
543 
544  DMatrix NK(m,m-1);
545  NK.setZero();
546 
547  for( run3 = 0; run3 < m-1; run3++ )
548  NK(run3,run3) = 1.0;
549 
550  for( run3 = 0; run3 < m-1; run3++ )
551  NK(m-1,run3) = -1.0;
552 
553 
554  DVector W(m);
555  for( run1 = 0; run1 < m; run1++ )
556  W(run1) = idx[run1];
557 
558  DVector PW = PHI_N*W;
559 
561 
562  Expression *Fnorm;
563  Fnorm = new Expression[m];
564 
565  for( run2 = 0; run2 < m; run2++ ){
566  Expression tmp3(1);
567  for( run3 = 0; run3 < m; run3++ ){
568  tmp3 = tmp3 + T(run2,run3)*( *arg[run3]- U(run3) );
569  }
570  Fnorm[run2] = tmp3;
571  }
572 
573  tmp.addMayerTerm( Fnorm[m-1] );
574  ocp_->setObjective( tmp );
575 
576  for( run1 = 0; run1 < m-1; run1++ ){
577 
578  Expression sum(1);
579 
580  for( run2 = 0; run2 < m; run2++ ){
581 
582  sum = sum + NK(run2,run1)*( Fnorm[run2] - PW(run2) );
583  }
584  ocp_->subjectTo( AT_END, sum <= 0.0 );
585  }
586  delete[] Fnorm;
587 
588  return SUCCESSFUL_RETURN;
589  }
590 
591 
592  // NORMAL BOUNDARY INTERSECTION:
593  // -----------------------------------------------------------------------------
594 
595  if( paretoGeneration == PFG_NORMAL_BOUNDARY_INTERSECTION ){
596 
597  DVector W(m);
598  for( run1 = 0; run1 < m; run1++ )
599  W(run1) = idx[run1];
600 
601  DMatrix P = getPayOffMatrix();
603  DVector V = P*W + U;
604  W.setAll( 1.0 );
605  DVector X = P*W;
606 
607  Expression lambda;
608 
609  lambda = ( *arg[m-1] - V(m-1) )/ X(m-1);
610 
611  tmp.addMayerTerm( *arg[m-1] );
612  ocp_->setObjective( tmp );
613 
614  for( run1 = 0; run1 < m-1; run1++ )
615  ocp->subjectTo( AT_END, *arg[run1] - lambda*X(run1) == V(run1) );
616 
617  return SUCCESSFUL_RETURN;
618  }
619  return SUCCESSFUL_RETURN;
620 }
621 
622 
624  VariablesGrid &xa_ ,
625  VariablesGrid &p_ ,
626  VariablesGrid &u_ ,
627  VariablesGrid &w_ ,
628  Expression **arg ){
629 
630  int run1;
631  Grid tmp_grid;
632  ocp->getGrid( tmp_grid );
633 
634  VariablesGrid *_xd = 0;
635  VariablesGrid *_xa = 0;
636  VariablesGrid *_p = 0;
637  VariablesGrid *_u = 0;
638  VariablesGrid *_w = 0;
639 
640  if( xd_.isEmpty() == BT_FALSE ) _xd = new VariablesGrid(xd_);
641  if( xa_.isEmpty() == BT_FALSE ) _xa = new VariablesGrid(xa_);
642  if( p_.isEmpty() == BT_FALSE ) _p = new VariablesGrid(p_ );
643  if( u_.isEmpty() == BT_FALSE ) _u = new VariablesGrid(u_ );
644  if( w_.isEmpty() == BT_FALSE ) _w = new VariablesGrid(w_ );
645 
646  Objective *obj;
647  for( run1 = 0; run1 < m; run1++ ){
648  obj = new Objective( tmp_grid );
649  obj->addMayerTerm(*arg[run1]);
650  OCPiterate xx( _xd, _xa, _p, _u, _w );
651  obj->evaluate( xx );
652  obj->getObjectiveValue( result(count,run1) );
653  delete obj;
654  }
655  count++;
656 
657  if( _xd != 0 ) delete _xd;
658  if( _xa != 0 ) delete _xa;
659  if( _p != 0 ) delete _p ;
660  if( _u != 0 ) delete _u ;
661  if( _w != 0 ) delete _w ;
662 
663  return SUCCESSFUL_RETURN;
664 }
665 
666 
667 
669 
670 // end of file.
const int defaultParetoFrontGeneration
virtual returnValue initializeNlpSolver(const OCPiterate &_userInit)
Data class for storing generic optimization variables.
Definition: ocp_iterate.hpp:57
USING_NAMESPACE_ACADO typedef TaylorVariable< Interval > T
VariablesGrid * x
DMatrix getWeights() const
void init(unsigned _nRows=0, unsigned _nCols=0)
Definition: matrix.hpp:135
returnValue getObjective(Objective &objective_) const
Definition: ocp.cpp:330
MultiObjectiveAlgorithm & operator=(const MultiObjectiveAlgorithm &arg)
const double INFTY
VariablesGrid * u
returnValue setStatus(BlockStatus _status)
DVector getNormalizationVector() const
BEGIN_NAMESPACE_ACADO const double EPS
User-interface to formulate and solve optimal control problems and static NLPs.
Generates weights for solving OCPs having multiple objectives.
Stores and evaluates the constraints of optimal control problems.
Definition: constraint.hpp:60
DMatrix getPayOffMatrix() const
Provides a time grid consisting of vector-valued optimization variables at each grid point...
Allows to pass back messages to the calling function.
void setAll(const T &_value)
Definition: matrix.hpp:141
returnValue addMayerTerm(const Expression &arg)
Allows to conveniently handle (one-dimensional) grids consisting of time points.
Definition: grid.hpp:58
BooleanType isEmpty() const
returnValue getParameters(VariablesGrid &u_) const
returnValue subjectTo(const DifferentialEquation &differentialEquation_)
Definition: ocp.cpp:153
double acadoGetTime()
#define CLOSE_NAMESPACE_ACADO
returnValue getConstraint(Constraint &constraint_) const
Definition: ocp.cpp:331
returnValue evaluateObjectives(VariablesGrid &xd_, VariablesGrid &xa_, VariablesGrid &p_, VariablesGrid &u_, VariablesGrid &w_, Expression **arg1)
virtual returnValue getObjectiveValue(double &objectiveValue)
Definition: objective.cpp:402
returnValue getAlgebraicStates(VariablesGrid &xa_) const
VariablesGrid * xa
returnValue setConstraint(const Constraint &constraint_)
Definition: ocp.cpp:339
Base class for all variables within the symbolic expressions family.
Definition: expression.hpp:56
const int defaultParetoFrontDiscretization
DVector getUtopiaVector() const
const int defaultParetoFrontHotstart
unsigned getDim() const
Definition: matrix.hpp:181
int getNumberOfSteps() const
returnValue getDisturbances(VariablesGrid &w_) const
Derived & setZero(Index size)
Data class for defining optimal control problems.
Definition: ocp.hpp:89
void hotstart(const real_t *const g, const real_t *const lb, const real_t *const ub, const real_t *const lbA, const real_t *const ubA, int nWSR, Options *options, int nOutputs, mxArray *plhs[])
unsigned getNumRows() const
Definition: matrix.hpp:185
#define L
virtual returnValue initializeObjective(Objective *F)
returnValue evaluate(const OCPiterate &x)
Definition: objective.cpp:254
returnValue getWeights(const int &m, const int &pnts, const DVector &weightsLB, const DVector &weightsUB, DMatrix &Weights, DVector &formers) const
returnValue getControls(VariablesGrid &p_) const
#define ASSERT(x)
DMatrix getUtopiaPlaneVectors() const
VariablesGrid * p
#define BT_TRUE
Definition: acado_types.hpp:47
virtual returnValue initializeNlpSolver(const OCPiterate &_userInit)
unsigned getNumCols() const
Definition: matrix.hpp:189
void setAll(const T &_value)
Definition: vector.hpp:160
virtual returnValue solveSingleObjective(const int &number)
returnValue formulateOCP(double *idx, OCP *ocp_, Expression **arg)
returnValue getGrid(Grid &grid_) const
Definition: ocp.cpp:329
DMatrix getNormalizedPayOffMatrix() const
VariablesGrid * w
#define BEGIN_NAMESPACE_ACADO
#define BT_FALSE
Definition: acado_types.hpp:49
User-interface to formulate and solve optimal control problems with multiple objectives.
returnValue addOption(OptionsName name, int value)
Definition: options.cpp:301
OptimizationAlgorithm & operator=(const OptimizationAlgorithm &arg)
returnValue setObjective(const Objective &objective_)
Definition: ocp.cpp:334
virtual returnValue solve()
Stores and evaluates the objective function of optimal control problems.
Definition: objective.hpp:123
returnValue getDifferentialStates(VariablesGrid &xd_) const
#define ACADOERROR(retval)


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