DOpE
ipopt_problem.h
Go to the documentation of this file.
1 
24 #ifndef _IPOPT_PROBLEM_H_
25 #define _IPOPT_PROBLEM_H_
26 
27 #ifdef WITH_IPOPT
28 #include "IpTNLP.hpp"
29 #endif
30 
31 #include "controlvector.h"
32 
33 #include <iostream>
34 
35 namespace DOpE
36 {
37 #ifdef WITH_IPOPT
38 
49  template <typename RPROBLEM, typename VECTOR>
50  class Ipopt_Problem : public Ipopt::TNLP
51  {
52  public:
53  Ipopt_Problem( int& ret_val,
54  RPROBLEM * OP, ControlVector<VECTOR>& q,
55  const ControlVector<VECTOR>* q_min,
56  const ControlVector<VECTOR>* q_max,
57  const ConstraintVector<VECTOR>& c);
58 
59  virtual ~Ipopt_Problem() {}
60 
64  virtual bool get_nlp_info(Ipopt::Index& n, Ipopt::Index& m, Ipopt::Index& nnz_jac_g,
65  Ipopt::Index& nnz_h_lag, IndexStyleEnum& index_style);
66 
68  virtual bool get_bounds_info(Ipopt::Index n, Ipopt::Number* x_l, Ipopt::Number* x_u,
69  Ipopt::Index m, Ipopt::Number* g_l, Ipopt::Number* g_u);
70 
72  virtual bool get_starting_point(Ipopt::Index n, bool init_x, Ipopt::Number* x,
73  bool init_z, Ipopt::Number* z_L, Ipopt::Number* z_U,
74  Ipopt::Index m, bool init_lambda,
75  Ipopt::Number* lambda);
76 
78  virtual bool eval_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number& obj_value);
79 
81  virtual bool eval_grad_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number* grad_f);
82 
84  virtual bool eval_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index m, Ipopt::Number* g);
85 
90  virtual bool eval_jac_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x,
91  Ipopt::Index m, Ipopt::Index nele_jac, Ipopt::Index* iRow, Ipopt::Index *jCol,
92  Ipopt::Number* values);
93 
98  virtual bool eval_h(Ipopt::Index n, const Ipopt::Number* x, bool new_x,
99  Ipopt::Number obj_factor, Ipopt::Index m, const Ipopt::Number* lambda,
100  bool new_lambda, Ipopt::Index nele_hess, Ipopt::Index* iRow,
101  Ipopt::Index* jCol, Ipopt::Number* values);
102 
104 
108  virtual void finalize_solution(Ipopt::SolverReturn status,
109  Ipopt::Index n, const Ipopt::Number* x, const Ipopt::Number* z_L, const Ipopt::Number* z_U,
110  Ipopt::Index m, const Ipopt::Number* g, const Ipopt::Number* lambda,
111  Ipopt::Number obj_value,
112  const Ipopt::IpoptData* ip_data,
113  Ipopt::IpoptCalculatedQuantities* ip_cq);
115  private:
116  int& _ret_val;
117  RPROBLEM* _P;
118  ControlVector<VECTOR> _q;
119  ControlVector<VECTOR>& _init; //Also the return value!
120  const ControlVector<VECTOR>* _q_min;
121  const ControlVector<VECTOR>* _q_max;
122  ConstraintVector<VECTOR> _c;
123  };
124  /***************************************************************************************/
125  /****************************************IMPLEMENTATION*********************************/
126  /***************************************************************************************/
127  template <typename RPROBLEM, typename VECTOR>
128  Ipopt_Problem<RPROBLEM,VECTOR>::Ipopt_Problem(int& ret_val,
129  RPROBLEM* OP, ControlVector<VECTOR>& q,
130  const ControlVector<VECTOR>* q_min,
131  const ControlVector<VECTOR>* q_max,
132  const ConstraintVector<VECTOR>& c)
133  : _ret_val(ret_val), _P(OP), _q(q), _init(q),
134  _q_min(q_min), _q_max(q_max), _c(c)
135  {
136  }
137 
138  template <typename RPROBLEM, typename VECTOR>
139  bool Ipopt_Problem<RPROBLEM,VECTOR>::get_nlp_info(Ipopt::Index& n, Ipopt::Index& m, Ipopt::Index& nnz_jac_g,
140  Ipopt::Index& nnz_h_lag, IndexStyleEnum& index_style)
141  {
142  n = _q.GetSpacialVector().size(); //n unknowns
143  m = _c.GetGlobalConstraints().size(); //Only Global constraints
144  nnz_jac_g = m*n; //Size of constraint jacobian
145  nnz_h_lag = 0; //Size of hessian! (n * n Don't compute!)
146  index_style = TNLP::C_STYLE; //C style indexing (0-based)
147  return true;
148  }
149 
150  template <typename RPROBLEM, typename VECTOR>
151  bool Ipopt_Problem<RPROBLEM,VECTOR>::get_bounds_info(Ipopt::Index n, Ipopt::Number* x_l, Ipopt::Number* x_u,
152  Ipopt::Index m, Ipopt::Number* g_l, Ipopt::Number* g_u)
153  {
154  assert(n == (int) _q.GetSpacialVector().size());
155  assert(m == (int) _c.GetGlobalConstraints().size());
156 
157  //Lower and upper bounds on q
158  const VECTOR& lb = _q_min->GetSpacialVector();
159  const VECTOR& ub = _q_max->GetSpacialVector();
160  for(int i = 0; i < n; i++)
161  {
162  x_l[i] = lb(i);
163  x_u[i] = ub(i);
164  }
165  //Global constraints are given such that feasible means <= 0
166  for(int i = 0; i < m; i++)
167  {
168  g_l[i] = -1.e+20;
169  g_u[i] = 0.;
170  }
171  return true;
172  }
173 
174  template <typename RPROBLEM, typename VECTOR>
175  bool Ipopt_Problem<RPROBLEM,VECTOR>::get_starting_point(Ipopt::Index n, bool init_x, Ipopt::Number* x,
176  bool init_z, Ipopt::Number* /*z_L*/, Ipopt::Number* /*z_U*/,
177  Ipopt::Index /*m*/, bool init_lambda,
178  Ipopt::Number* /*lambda*/)
179  {
180  assert(init_x == true);
181  assert(init_z == false);
182  assert(init_lambda == false);
183  const VECTOR& in = _init.GetSpacialVector();
184 
185  for(int i = 0; i < n; i++)
186  {
187  x[i] = in(i);
188  }
189  return true;
190  }
191 
192  template <typename RPROBLEM, typename VECTOR>
193  bool Ipopt_Problem<RPROBLEM,VECTOR>::eval_f(Ipopt::Index n, const Ipopt::Number* x, bool /*new_x*/, Ipopt::Number& obj_value)
194  {
195  VECTOR& qval = _q.GetSpacialVector();
196  for(int i = 0; i < n; i++)
197  {
198  qval(i) = x[i];
199  }
200  try
201  {
202  obj_value = _P->ComputeReducedCostFunctional(_q);
203  }
204  catch(DOpEException& e)
205  {
206  _P->GetExceptionHandler()->HandleCriticalException(e,"IPOPT_RPROBLEM::eval_f");
207  }
208  return true;
209  }
210 
211  template <typename RPROBLEM, typename VECTOR>
212  bool Ipopt_Problem<RPROBLEM,VECTOR>::eval_grad_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number* grad_f)
213  {
214  ControlVector<VECTOR> gradient(_q);
215  ControlVector<VECTOR> gradient_transposed(_q);
216  if( new_x )
217  {
218  //Need to calculate J!
219  VECTOR& qval = _q.GetSpacialVector();
220  for(int i = 0; i < n; i++)
221  {
222  qval(i) = x[i];
223  }
224  try
225  {
226  _P->ComputeReducedCostFunctional(_q);
227  }
228  catch(DOpEException& e)
229  {
230  _P->GetExceptionHandler()->HandleCriticalException(e,"IPOPT_PROBLEM::eval_grad_f");
231  }
232  }
233  //Compute Functional Gradient
234  try
235  {
236  _P->ComputeReducedGradient(_q,gradient,gradient_transposed);
237  }
238  catch(DOpEException& e)
239  {
240  _P->GetExceptionHandler()->HandleCriticalException(e,"IPOPT_PROBLEM::eval_grad_f");
241  }
242  const VECTOR& ref_g = gradient_transposed.GetSpacialVector();
243  for(int i=0; i < n; i++)
244  {
245  grad_f[i] = ref_g(i);
246  }
247 
248  return true;
249  }
250 
251  template <typename RPROBLEM, typename VECTOR>
252  bool Ipopt_Problem<RPROBLEM,VECTOR>::eval_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index /*m*/, Ipopt::Number* g)
253  {
254  if( new_x )
255  {
256  //Need to calculate J!
257  VECTOR& qval = _q.GetSpacialVector();
258  for(int i = 0; i < n; i++)
259  {
260  qval(i) = x[i];
261  }
262  try
263  {
264  _P->ComputeReducedCostFunctional(_q);
265  }
266  catch(DOpEException& e)
267  {
268  _P->GetExceptionHandler()->HandleCriticalException(e,"IPOPT_PROBLEM::eval_g");
269  }
270  }
271  //Calculate constraints
272  try
273  {
274  _P->ComputeReducedConstraints(_q,_c);
275  }
276  catch(DOpEException& e)
277  {
278  _P->GetExceptionHandler()->HandleCriticalException(e,"IPOPT_PROBLEM::eval_g");
279  }
280  const dealii::Vector<double>& gc = _c.GetGlobalConstraints();
281  for(unsigned int i=0; i < gc.size(); i++)
282  {
283  g[i] = gc(i);
284  }
285  return true;
286  }
287 
288  template <typename RPROBLEM, typename VECTOR>
289  bool Ipopt_Problem<RPROBLEM,VECTOR>::eval_jac_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x,
290  Ipopt::Index m, Ipopt::Index nele_jac, Ipopt::Index* iRow, Ipopt::Index *jCol,
291  Ipopt::Number* values)
292  {
293  if(values != NULL)
294  {
295  //Compute Jacobian of Constraints
296  ControlVector<VECTOR> gradient(_q);
297  ControlVector<VECTOR> gradient_transposed(_q);
298  if( new_x )
299  {
300  //Need to calculate J!
301  VECTOR& qval = _q.GetSpacialVector();
302  for(int i = 0; i < n; i++)
303  {
304  qval(i) = x[i];
305  }
306  try
307  {
308  _P->ComputeReducedCostFunctional(_q);
309  }
310  catch(DOpEException& e)
311  {
312  _P->GetExceptionHandler()->HandleCriticalException(e,"IPOPT_PROBLEM::eval_grad_g");
313  }
314  }
315  //Compute Constraint Gradients
316  for(int j=0; j < m; j++)
317  {
318  try
319  {
320  _P->ComputeReducedGradientOfGlobalConstraints(j,_q,_c,
321  gradient,gradient_transposed);
322  }
323  catch(DOpEException& e)
324  {
325  _P->GetExceptionHandler()->HandleCriticalException(e,"IPOPT_PROBLEM::eval_grad_g");
326  }
327  const VECTOR& ref_g = gradient_transposed.GetSpacialVector();
328  for(int i=0; i < n; i++)
329  {
330  values[n*j+i] = ref_g(i);
331  }
332  }
333  }
334  else
335  {
336  assert(nele_jac == n*m);
337  // return the structure of the Jacobian (here a dense one)
338  for(int i = 0; i < n; i++)
339  {
340  for(int j = 0; j < m; j++)
341  {
342  iRow[n*j+i] = j;
343  jCol[n*j+i] = i;
344  }
345  }
346  }
347 
348  return true;
349  }
350 
351  template <typename RPROBLEM, typename VECTOR>
352  bool Ipopt_Problem<RPROBLEM,VECTOR>::eval_h(Ipopt::Index /*n*/, const Ipopt::Number* /*x*/, bool /*new_x*/,
353  Ipopt::Number /*obj_factor*/, Ipopt::Index /*m*/, const Ipopt::Number* /*lambda*/,
354  bool /*new_lambda*/, Ipopt::Index /*nele_hess*/, Ipopt::Index* /*iRow*/,
355  Ipopt::Index* /*jCol*/, Ipopt::Number* /*values*/)
356  {
357  //Donot calculate the hessian!
358  return false;
359  }
360 
361 
362  template <typename RPROBLEM, typename VECTOR>
363  void Ipopt_Problem<RPROBLEM,VECTOR>::finalize_solution(Ipopt::SolverReturn status,
364  Ipopt::Index n, const Ipopt::Number* x,
365  const Ipopt::Number* /*z_L*/, const Ipopt::Number* /*z_U*/,
366  Ipopt::Index /*m*/, const Ipopt::Number* /*g*/,
367  const Ipopt::Number* /*lambda*/,
368  Ipopt::Number /*obj_value*/,
369  const Ipopt::IpoptData* /*ip_data*/,
370  Ipopt::IpoptCalculatedQuantities* /*ip_cq*/)
371  {
372  //Check the result, _q should be x and then be copied to _init
373 
374  VECTOR& ret = _init.GetSpacialVector();
375  for(int i = 0; i < n; i++)
376  {
377  ret(i) = x[i];
378  }
379  if(status == Ipopt::SUCCESS)
380  _ret_val = 1;
381  else
382  _ret_val = 0;
383  }
384 
385 #endif //Endof WITH_IPOPT
386 
387 
388 
389 } //Endof Namespace DOpE
390 
391 #endif
Generated on Fri Oct 24 2014 17:56:40 for DOpE by   doxygen 1.8.6