FFSM++  1.1.0 French Forest Sector Model ++
Ipopt_nlp_problem_debugtest.cpp
Go to the documentation of this file.
2
3 #include <cassert>
4 #include <iostream>
5
6 using namespace Ipopt;
7
8 // constructor
10 {}
11
12 //destructor
14 {}
15
16 // returns the size of the problem
17 bool Ipopt_nlp_problem_debugtest::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
18  Index& nnz_h_lag, IndexStyleEnum& index_style)
19 {
20  // The problem described in Ipopt_nlp_problem_debugtest.hpp has 4 variables, x[0] through x[3]
21  n = 4;
22
23  // one equality constraint and one inequality constraint
24  m = 2;
25
26  // in this example the jacobian is dense and contains 8 nonzeros
27  nnz_jac_g = 8;
28
29  // the hessian is also dense and has 16 total nonzeros, but we
30  // only need the lower left corner (since it is symmetric)
31  nnz_h_lag = 10;
32
33  // use the C style indexing (0-based)
34  index_style = TNLP::C_STYLE;
35
36  return true;
37 }
38
39 // returns the variable bounds
40 bool Ipopt_nlp_problem_debugtest::get_bounds_info(Index n, Number* x_l, Number* x_u,
41  Index m, Number* g_l, Number* g_u)
42 {
43  // here, the n and m we gave IPOPT in get_nlp_info are passed back to us.
44  // If desired, we could assert to make sure they are what we think they are.
45  assert(n == 4);
46  assert(m == 2);
47
48  // the variables have lower bounds of 1
49  for (Index i=0; i<4; i++) {
50  x_l[i] = 1.0;
51  }
52
53  // the variables have upper bounds of 5
54  for (Index i=0; i<4; i++) {
55  x_u[i] = 5.0;
56  }
57
58  // the first constraint g1 has a lower bound of 25
59  g_l[0] = 25;
60  // the first constraint g1 has NO upper bound, here we set it to 2e19.
61  // Ipopt interprets any number greater than nlp_upper_bound_inf as
62  // infinity. The default value of nlp_upper_bound_inf and nlp_lower_bound_inf
63  // is 1e19 and can be changed through ipopt options.
64  g_u[0] = 2e19;
65
66  // the second constraint g2 is an equality constraint, so we set the
67  // upper and lower bound to the same value
68  g_l[1] = g_u[1] = 40.0;
69
70  return true;
71 }
72
73 // returns the initial point for the problem
74 bool Ipopt_nlp_problem_debugtest::get_starting_point(Index n, bool init_x, Number* x,
75  bool init_z, Number* z_L, Number* z_U,
76  Index m, bool init_lambda,
77  Number* lambda)
78 {
79  // Here, we assume we only have starting values for x, if you code
80  // your own NLP, you can provide starting values for the dual variables
81  // if you wish
82  assert(init_x == true);
83  assert(init_z == false);
84  assert(init_lambda == false);
85
86  // initialize to the given starting point
87  x[0] = 1.0;
88  x[1] = 5.0;
89  x[2] = 5.0;
90  x[3] = 1.0;
91
92  return true;
93 }
94
95 // returns the value of the objective function
96 bool Ipopt_nlp_problem_debugtest::eval_f(Index n, const Number* x, bool new_x, Number& obj_value)
97 {
98  assert(n == 4);
99
100  obj_value = x[0] * x[3] * (x[0] + x[1] + x[2]) + x[2];
101
102  return true;
103 }
104
107 {
108  assert(n == 4);
109
110  grad_f[0] = x[0] * x[3] + x[3] * (x[0] + x[1] + x[2]);
111  grad_f[1] = x[0] * x[3];
112  grad_f[2] = x[0] * x[3] + 1;
113  grad_f[3] = x[0] * (x[0] + x[1] + x[2]);
114
115  return true;
116 }
117
118 // return the value of the constraints: g(x)
119 bool Ipopt_nlp_problem_debugtest::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
120 {
121  assert(n == 4);
122  assert(m == 2);
123
124  g[0] = x[0] * x[1] * x[2] * x[3];
125  g[1] = x[0]*x[0] + x[1]*x[1] + x[2]*x[2] + x[3]*x[3];
126
127  return true;
128 }
129
130 // return the structure or values of the jacobian
131 bool Ipopt_nlp_problem_debugtest::eval_jac_g(Index n, const Number* x, bool new_x,
132  Index m, Index nele_jac, Index* iRow, Index *jCol,
133  Number* values)
134 {
135  if (values == NULL) {
136  // return the structure of the jacobian
137
138  // this particular jacobian is dense
139  iRow[0] = 0;
140  jCol[0] = 0;
141  iRow[1] = 0;
142  jCol[1] = 1;
143  iRow[2] = 0;
144  jCol[2] = 2;
145  iRow[3] = 0;
146  jCol[3] = 3;
147  iRow[4] = 1;
148  jCol[4] = 0;
149  iRow[5] = 1;
150  jCol[5] = 1;
151  iRow[6] = 1;
152  jCol[6] = 2;
153  iRow[7] = 1;
154  jCol[7] = 3;
155  }
156  else {
157  // return the values of the jacobian of the constraints
158
159  values[0] = x[1]*x[2]*x[3]; // 0,0
160  values[1] = x[0]*x[2]*x[3]; // 0,1
161  values[2] = x[0]*x[1]*x[3]; // 0,2
162  values[3] = x[0]*x[1]*x[2]; // 0,3
163
164  values[4] = 2*x[0]; // 1,0
165  values[5] = 2*x[1]; // 1,1
166  values[6] = 2*x[2]; // 1,2
167  values[7] = 2*x[3]; // 1,3
168  }
169
170  return true;
171 }
172
173
174 //return the structure or values of the hessian
175 bool Ipopt_nlp_problem_debugtest::eval_h(Index n, const Number* x, bool new_x,
176  Number obj_factor, Index m, const Number* lambda,
177  bool new_lambda, Index nele_hess, Index* iRow,
178  Index* jCol, Number* values)
179 {
180  if (values == NULL) {
181  // return the structure. This is a symmetric matrix, fill the lower left
182  // triangle only.
183
184  // the hessian for this problem is actually dense
185  Index idx=0;
186  for (Index row = 0; row < 4; row++) {
187  for (Index col = 0; col <= row; col++) {
188  iRow[idx] = row;
189  jCol[idx] = col;
190  idx++;
191  }
192  }
193
194  assert(idx == nele_hess);
195  }
196  else {
197  // return the values. This is a symmetric matrix, fill the lower left
198  // triangle only
199
200  // fill the objective portion
201  values[0] = obj_factor * (2*x[3]); // 0,0
202
203  values[1] = obj_factor * (x[3]); // 1,0
204  values[2] = 0.; // 1,1
205
206  values[3] = obj_factor * (x[3]); // 2,0
207  values[4] = 0.; // 2,1
208  values[5] = 0.; // 2,2
209
210  values[6] = obj_factor * (2*x[0] + x[1] + x[2]); // 3,0
211  values[7] = obj_factor * (x[0]); // 3,1
212  values[8] = obj_factor * (x[0]); // 3,2
213  values[9] = 0.; // 3,3
214
215
216  // add the portion for the first constraint
217  values[1] += lambda[0] * (x[2] * x[3]); // 1,0
218
219  values[3] += lambda[0] * (x[1] * x[3]); // 2,0
220  values[4] += lambda[0] * (x[0] * x[3]); // 2,1
221
222  values[6] += lambda[0] * (x[1] * x[2]); // 3,0
223  values[7] += lambda[0] * (x[0] * x[2]); // 3,1
224  values[8] += lambda[0] * (x[0] * x[1]); // 3,2
225
226  // add the portion for the second constraint
227  values[0] += lambda[1] * 2; // 0,0
228
229  values[2] += lambda[1] * 2; // 1,1
230
231  values[5] += lambda[1] * 2; // 2,2
232
233  values[9] += lambda[1] * 2; // 3,3
234  }
235
236  return true;
237 }
238
239
240
242  Index n, const Number* x, const Number* z_L, const Number* z_U,
243  Index m, const Number* g, const Number* lambda,
244  Number obj_value,
245  const IpoptData* ip_data,
246  IpoptCalculatedQuantities* ip_cq)
247 {
248  // here is where we would store the solution to variables, or write to a file, etc
249  // so we could use the solution.
250
251  // For this example, we write the solution to the console
252  std::cout << std::endl << std::endl << "Solution of the primal variables, x" << std::endl;
253  for (Index i=0; i<n; i++) {
254  std::cout << "x[" << i << "] = " << x[i] << std::endl;
255  }
256
257  std::cout << std::endl << std::endl << "Solution of the bound multipliers, z_L and z_U" << std::endl;
258  for (Index i=0; i<n; i++) {
259  std::cout << "z_L[" << i << "] = " << z_L[i] << std::endl;
260  }
261  for (Index i=0; i<n; i++) {
262  std::cout << "z_U[" << i << "] = " << z_U[i] << std::endl;
263  }
264
265  std::cout << std::endl << std::endl << "Objective value" << std::endl;
266  std::cout << "f(x*) = " << obj_value << std::endl;
267
268  std::cout << std::endl << "Final value of the constraints:" << std::endl;
269  for (Index i=0; i<m ;i++) {
270  std::cout << "g(" << i << ") = " << g[i] << std::endl;
271  }
272 }
virtual bool eval_h(Index n, const Number *x, bool new_x, Number obj_factor, Index m, const Number *lambda, bool new_lambda, Index nele_hess, Index *iRow, Index *jCol, Number *values)
virtual bool eval_f(Index n, const Number *x, bool new_x, Number &obj_value)
virtual bool get_bounds_info(Index n, Number *x_l, Number *x_u, Index m, Number *g_l, Number *g_u)
virtual bool get_nlp_info(Index &n, Index &m, Index &nnz_jac_g, Index &nnz_h_lag, IndexStyleEnum &index_style)
virtual void finalize_solution(SolverReturn status, Index n, const Number *x, const Number *z_L, const Number *z_U, Index m, const Number *g, const Number *lambda, Number obj_value, const IpoptData *ip_data, IpoptCalculatedQuantities *ip_cq)
virtual bool eval_jac_g(Index n, const Number *x, bool new_x, Index m, Index nele_jac, Index *iRow, Index *jCol, Number *values)
virtual bool get_starting_point(Index n, bool init_x, Number *x, bool init_z, Number *z_L, Number *z_U, Index m, bool init_lambda, Number *lambda)