source:branches/dev/Examples/hs071_cpp/hs071_nlp.cpp@461

Last change on this file since 461 was 461, checked in by claird, 9 years ago

Fixed error in the example and the documentation, wrong starting point.

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