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