/[ascend]/trunk/solvers/ipopt/asc_ipopt.c
ViewVC logotype

Contents of /trunk/solvers/ipopt/asc_ipopt.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2351 - (show annotations) (download) (as text)
Thu Jan 6 02:02:41 2011 UTC (11 years, 6 months ago) by jpye
File MIME type: text/x-csrc
File size: 46287 byte(s)
Resolved memory leak with test/test solver_ipopt.formula.
Still some issues to overcome with regard to inactive relations, probably.
Some issue still occurring when attempting to run whole solver_ipopt suite in one go.
1 /* ASCEND modelling environment
2 Copyright (C) 2007-2010 Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330,
17 Boston, MA 02111-1307, USA.
18 *//**
19 @file
20 Connection of the IPOPT optimisation solver into ASCEND.
21
22 THIS IS STILL VERY MUCH UNDER DEVELOPMENT AND INCOMPLETE.
23
24 The IPOPT solver is documented at http://projects.coin-or.org/Ipopt/
25 *//*
26 ASCEND wrapper for IPOPT originally by John Pye, Jun 2007 onwards. Further
27 contributions from Mahesh Narayanamurthi 2009.
28 */
29
30 #include <ascend/utilities/config.h>
31
32 #ifndef ASC_WITH_IPOPT
33 # error "ASC_WITH_IPOPT must be defined in order to build this."
34 #endif
35
36 #include <math.h>
37
38 #include <ascend/solver/solver.h>
39
40 #include <ascend/system/calc.h>
41 #include <ascend/system/relman.h>
42 #include <ascend/system/slv_stdcalls.h>
43 #include <ascend/system/block.h>
44
45 #include <ascend/general/platform.h>
46 #include <ascend/general/panic.h>
47 #include <ascend/general/ascMalloc.h>
48 #include <ascend/utilities/ascDynaLoad.h>
49 #include <ascend/general/mem.h>
50 #include <ascend/utilities/ascEnvVar.h>
51 #include <ascend/general/tm_time.h>
52 #include <ascend/general/env.h>
53
54 #include <ascend/general/ltmatrix.h>
55
56 #include <coin/IpStdCInterface.h>
57
58 ASC_DLLSPEC SolverRegisterFn ipopt_register;
59
60 /*------------------------------------------------------------------------------
61 DATA STRUCTURES AND FORWARD DEFS
62 */
63
64 /**
65 Documentation of solver options for IPOPT is at
66 http://www.coin-or.org/Ipopt/documentation/node1.html
67 */
68 enum{
69 /** ASCEND OPTIONS */
70 ASCEND_PARAM_SAFEEVAL
71 /** OUTPUT OPTIONS */
72 ,IPOPT_PARAM_PRINT_LEVEL
73 ,IPOPT_PARAM_PRINT_USER_OPTIONS
74 /** TERMINATION OPTIONS */
75 ,IPOPT_PARAM_TOL
76 ,IPOPT_PARAM_MAX_ITER
77 ,IPOPT_PARAM_MAX_CPU_TIME
78 ,IPOPT_PARAM_DIVERGING_ITERATES_TOL
79 ,IPOPT_PARAM_CONSTR_VIOL_TOL
80 ,IPOPT_PARAM_DUAL_INFEASIBILITY_TOL
81 ,IPOPT_PARAM_ACCEPTABLE_TOL
82 ,IPOPT_PARAM_ACCEPTABLE_ITER
83 /** LINEAR SOLVER OPTIONS */
84 ,IPOPT_PARAM_LINEAR_SOLVER
85 /** BARRIER PARAMETER OPTIONS */
86 ,IPOPT_PARAM_MU_STRATEGY
87 /** DERIVATIVE TEST OPTIONS */
88 ,IPOPT_PARAM_DERIVATIVE_TEST
89 /** QUASI-NEWTON OPTIONS */
90 ,IPOPT_PARAM_HESS_APPROX
91 /** OPTIONS COUNT */
92 ,IPOPT_PARAMS
93 };
94
95 #define SYS(s) ((IpoptSystem *)(s))
96
97 struct IpoptSystemStruct{
98
99 /*
100 Problem definition
101 */
102 slv_system_t slv; /* slv_system_t back-link */
103 struct rel_relation *obj; /* Objective function: NULL = none */
104 struct rel_relation *old_obj;/* Objective function: NULL = none */
105 struct var_variable **vlist; /* Variable list (NULL terminated) */
106 struct rel_relation **rlist; /* Relation list (NULL terminated) */
107
108 var_filter_t vfilt;
109 rel_filter_t rfilt;
110
111 /*
112 Solver information
113 */
114 int32 presolved; /* ? Has the system been presolved */
115 int32 resolve; /* ? Has the system been resolved */
116 slv_parameters_t p; /* Parameters */
117 slv_status_t s; /* Status (as of iteration end) */
118
119 int32 cap; /* Order of matrix/vectors */
120 int32 rank; /* Symbolic rank of problem */
121 int32 vused; /* Free and incident variables */
122 int32 vtot; /* length of varlist */
123 int32 rused; /* Included relations */
124 int32 rtot; /* length of rellist */
125 double clock; /* CPU time */
126
127 int32 calc_ok;
128 double obj_val;
129
130 #if 0
131 void *parm_array[IPOPT_PARAMS];
132 #endif
133 struct slv_parameter pa[IPOPT_PARAMS];
134
135 /*
136 IPOPT DATA
137 */
138 Index n; /* number of variables */
139 Index m; /* number of constraints (excl the 'objective relation')*/
140
141 Index nnzJ; /* number of non zeros in the jacobian of the constraints */
142 Index nnzH; /* number of non-zeros in the hessian of the objective */
143
144 #if 0
145 Number* x_L; /* lower bounds on x */
146 Number* x_U; /* upper bounds on x */
147 Number* g_L; /* lower bounds on g */
148 Number* g_U; /* upper bounds on g */
149 #endif
150
151 IpoptProblem nlp; /* IpoptProblem */
152
153 enum ApplicationReturnStatus status; /* Solve return code */
154 #if 0
155 Number* x; /* starting point and solution vector */
156 Number* mult_x_L; /* lower bound multipliers at the solution */
157 Number* mult_x_U; /* upper bound multipliers at the solution */
158 #endif
159 Index i; /* generic counter */
160 };
161
162 typedef struct IpoptSystemStruct IpoptSystem;
163
164 static int ipopt_get_default_parameters(slv_system_t server, SlvClientToken asys, slv_parameters_t *parameters);
165
166 static void ipopt_iteration_begins(IpoptSystem *sys);
167 static void ipopt_iteration_ends(IpoptSystem *sys);
168
169 /*------------------------------------------------------------------------------
170 SYSTEM SETUP/DESTROY, STATUS AND SOLVER ELIGIBILITY
171 */
172
173 static SlvClientToken ipopt_create(slv_system_t server, int32 *statusindex){
174 IpoptSystem *sys;
175
176 sys = ASC_NEW_CLEAR(IpoptSystem);
177 if(sys==NULL){
178 *statusindex = 1;
179 return sys;
180 }
181
182 sys->p.parms = sys->pa;
183 sys->p.dynamic_parms = 0;
184 ipopt_get_default_parameters(server,(SlvClientToken)sys,&(sys->p));
185 sys->p.whose = (*statusindex);
186
187 sys->presolved = 0;
188 sys->resolve = 0;
189
190 sys->n = -1;
191 sys->m = -1;
192
193 sys->s.ok = TRUE;
194 sys->s.calc_ok = TRUE;
195 sys->s.costsize = 0;
196 sys->s.cost = NULL; /*redundant, but sanity preserving */
197 sys->s.block.number_of = 1;
198 sys->s.block.current_block = 0;
199 sys->s.block.current_reordered_block = 0;
200 sys->s.block.current_size = 0;
201 sys->s.block.previous_total_size = 0;
202 sys->s.block.iteration = 0;
203 sys->s.block.funcs = 0;
204 sys->s.block.jacs = 0;
205 sys->s.block.cpu_elapsed = 0;
206 sys->s.block.functime = 0;
207 sys->s.block.jactime = 0;
208 sys->s.block.residual = 0;
209
210 sys->rfilt.matchbits = (REL_INCLUDED | REL_ACTIVE);
211 sys->rfilt.matchvalue = (REL_INCLUDED | REL_ACTIVE);
212 sys->vfilt.matchbits = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR | VAR_FIXED);
213 sys->vfilt.matchvalue = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR);
214
215 sys->vlist = slv_get_solvers_var_list(server);
216 sys->rlist = slv_get_solvers_rel_list(server);
217
218 sys->rtot = slv_get_num_solvers_rels(server);
219 sys->vtot = slv_get_num_solvers_vars(server);
220
221 sys->obj = slv_get_obj_relation(server);
222
223 sys->slv = server;
224
225 /*char *tmp = rel_make_name(sys->slv,sys->obj);
226 //CONSOLE_DEBUG("Objective relation is '%s'",tmp);
227 ASC_FREE(tmp);*/
228
229 //CONSOLE_DEBUG("There are %d constraint relations.", sys->rtot);
230
231 if(sys->vlist == NULL) {
232 ASC_FREE(sys);
233 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IPOPT called with no variables.");
234 *statusindex = -2;
235 return NULL;
236 }
237
238 if(sys->rlist == NULL && sys->obj == NULL) {
239 ASC_FREE(sys);
240 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IPOPT called with no relations or objective.");
241 *statusindex = -1;
242 return NULL;
243 }
244
245
246 /* do nothing with the objective list, pars, bounds, extrels, etc etc */
247
248 slv_check_var_initialization(server);
249 *statusindex = 0;
250 return((SlvClientToken)sys);
251 }
252
253 static int32 ipopt_destroy(slv_system_t server, SlvClientToken asys){
254 IpoptSystem *sys;
255 UNUSED_PARAMETER(server);
256 sys = SYS(asys);
257 slv_destroy_parms(&(sys->p));
258 if(sys->s.cost) ascfree(sys->s.cost);
259 ASC_FREE(sys);
260 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"ipopt_destroy still needs debugging");
261 return 0;
262 }
263
264
265 static int ipopt_get_status(slv_system_t server, SlvClientToken asys
266 ,slv_status_t *status
267 ){
268 IpoptSystem *sys;
269 (void)server; /* stop gcc whine about unused parameter */
270
271 sys = SYS(asys);
272 //if (check_system(sys)) return 1;
273 mem_copy_cast(&(sys->s),status,sizeof(slv_status_t));
274 return 0;
275 }
276
277 /**
278 Update the solver status. FIXME can't we get rid of this silly function
279 somehot?
280 */
281 static void update_status(IpoptSystem *sys){
282 boolean unsuccessful;
283
284 sys->s.time_limit_exceeded = FALSE; /* can't do this one with IPOPT */
285 sys->s.iteration_limit_exceeded = FALSE; /* IPOPT handles this one internally */
286
287 unsuccessful = sys->s.diverged || sys->s.inconsistent ||
288 sys->s.iteration_limit_exceeded || sys->s.time_limit_exceeded;
289
290 sys->s.ready_to_solve = !unsuccessful && !sys->s.converged;
291 sys->s.ok = !unsuccessful && sys->s.calc_ok && !sys->s.struct_singular;
292 }
293
294 static int32 ipopt_eligible_solver(slv_system_t server){
295 struct rel_relation **rp;
296 struct var_variable **vp;
297 rel_filter_t rfilt;
298 var_filter_t vfilt;
299
300 rfilt.matchbits = (REL_CONDITIONAL | REL_INWHEN);
301 rfilt.matchvalue = (REL_CONDITIONAL | REL_INWHEN);
302
303 vfilt.matchbits = (VAR_BINARY);
304 vfilt.matchvalue = (VAR_BINARY);
305
306 /// @todo check that there is a MAXIMIZE or MINIMIZE statement
307 if (slv_get_obj_relation(server) == NULL)
308 ERROR_REPORTER_HERE(ASC_USER_ERROR,"No objective function found");
309
310 /// @todo check that there are no WHENs or CONDITIONALs
311 for( rp=slv_get_solvers_rel_list(server); *rp != NULL ; ++rp ) {
312 if(rel_apply_filter(*rp,&rfilt)){
313 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"WHEN and CONDITIONAL Statements are not supported.");
314 return(FALSE);
315 }
316 }
317
318 /// @todo check that there are no discrete-valued variables
319 for( vp=slv_get_solvers_var_list(server); *vp != NULL ; ++vp ) {
320 if(var_apply_filter(*vp,&vfilt)){
321 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Discrete Variables not supported.");
322 return(FALSE);
323 }
324 }
325
326 /// @todo check anything else?
327
328 return 1;
329 }
330
331 /*------------------------------------------------------------------------------
332 SOLVER PARAMETERS
333 */
334
335 static
336 int32 ipopt_get_default_parameters(slv_system_t server, SlvClientToken asys
337 ,slv_parameters_t *parameters
338 ){
339 IpoptSystem *sys = NULL;
340 struct slv_parameter *new_parms = NULL;
341 int32 make_macros = 0;
342
343 if(server != NULL && asys != NULL) {
344 sys = SYS(asys);
345 make_macros = 1;
346 }
347
348 if(parameters->parms == NULL) {
349 new_parms = ASC_NEW_ARRAY_OR_NULL(struct slv_parameter,IPOPT_PARAMS);
350 if(new_parms == NULL) {
351 return -1;
352 }
353 parameters->parms = new_parms;
354 parameters->dynamic_parms = 1;
355 }
356
357 parameters->num_parms = 0;
358
359 /** ASCEND Options */
360
361 slv_param_bool(parameters,ASCEND_PARAM_SAFEEVAL
362 ,(SlvParameterInitBool){{"safeeval"
363 ,"Use safe evaluation?",1
364 ,"Use 'safe' function evaluation routines (TRUE) or allow ASCEND to "
365 "throw SIGFPE errors which will then halt integration (FALSE)."
366 }, FALSE}
367 );
368
369
370 /** Output Options */
371
372 slv_param_int(parameters,IPOPT_PARAM_PRINT_LEVEL
373 ,(SlvParameterInitInt){{"print_level"
374 ,"Output verbosity level",2
375 ,"Sets the default verbosity level for console output."
376 " The larger this value the more detailed is the output."
377 " Default value is 5 and range is 0 to 12."
378 }, 5, 0, 12}
379 );
380
381
382 slv_param_char(parameters,IPOPT_PARAM_PRINT_USER_OPTIONS
383 ,(SlvParameterInitChar){{"print_user_options"
384 ,"Print all options set by the user.",2
385 ,"If selected, the algorithm will print the list of all options"
386 " set by the user including their values and whether they have"
387 " been used. In some cases this information might be incorrect,"
388 " due to the internal program flow. The default value for this "
389 " string option is 'no'. "
390 }, "yes"}, (char *[]){
391 "no","yes",NULL
392 }
393 );
394
395
396
397 /** Termination Options */
398
399 slv_param_int(parameters,IPOPT_PARAM_MAX_ITER
400 ,(SlvParameterInitInt){{"max_iter"
401 ,"Maximum number of iterations",3
402 ,"The algorithm terminates with an error message if the number of iterations exceeded this number."
403 }, 3000, 0, 100000000}
404 );
405
406 slv_param_real(parameters,IPOPT_PARAM_TOL
407 ,(SlvParameterInitReal){{"tol"
408 ,"Desired convergence tolerance (relative)",3
409 ,"Determines the convergence tolerance for the algorithm. The algorithm"
410 " terminates successfully, if the (scaled) NLP error becomes smaller"
411 " than this value, and if the (absolute) criteria according to "
412 " 'dual_inf_tol', 'primal_inf_tol', and 'cmpl_inf_tol' are met. (This"
413 " is epsilon_tol in Eqn. (6) in implementation paper). See also "
414 " 'acceptable_tol' as a second termination criterion. Note, some other"
415 " algorithmic features also use this quantity to determine thresholds"
416 " etc."
417 }, 1.e-8, 0, 1.e20}
418 );
419
420 slv_param_real(parameters,IPOPT_PARAM_MAX_CPU_TIME
421 ,(SlvParameterInitReal){{"max_cpu_time"
422 ,"Maximum CPU time allowed per problem (seconds)",3
423 ,"The algorithm terminates with an error message if the CPU time exceeds this value."
424 }, 1.e6, 0, 1.e7}
425 );
426
427 slv_param_real(parameters,IPOPT_PARAM_DIVERGING_ITERATES_TOL
428 ,(SlvParameterInitReal){{"diverging_iterates_tol"
429 ,"Threshold for maximal value of primal iterates",3
430 ,"If any component of the primal iterates exceeded this value"
431 " (in absolute terms), the optimization is aborted with the "
432 "exit message that the iterates seem to be diverging"
433 }, 1.e20, 0, 1.e50}
434 );
435
436
437 slv_param_real(parameters,IPOPT_PARAM_DUAL_INFEASIBILITY_TOL
438 ,(SlvParameterInitReal){{"dual_inf_tol"
439 ,"Desired threshold for the dual infeasibility.",3
440 ,"Absolute tolerance on the dual infeasibility. Successful "
441 "termination requires that the max-norm of the (unscaled) "
442 "dual infeasibility is less than this threshold. The valid "
443 "range for this real option is 0 < dual_inf_tol < +inf and"
444 " its default value is 1."
445 }, 1, 0, 1.e50}
446 );
447
448
449 slv_param_real(parameters,IPOPT_PARAM_CONSTR_VIOL_TOL
450 ,(SlvParameterInitReal){{"constr_viol_tol"
451 ,"Desired threshold for the constraint violation.",3
452 ,"Absolute tolerance on the constraint violation. Successful"
453 "termination requires that the max-norm of the (unscaled) "
454 " constraint violation is less than this threshold. The valid"
455 " range for this real option is 0 < constr_viol_tol < +inf and"
456 " its default value is 0.0001"
457 }, 1e-4, 0, 1.e50}
458 );
459
460 slv_param_real(parameters,IPOPT_PARAM_ACCEPTABLE_TOL
461 ,(SlvParameterInitReal){{"acceptable_tol"
462 ,"Acceptable convergence tolerance (relative).",3
463 ,"Determines which (scaled) overall optimality error is"
464 " considered to be 'acceptable.' There are two levels of"
465 " termination criteria. If the usual 'desired' tolerances"
466 " (see tol, dual_inf_tol etc) are satisfied at an iteration,"
467 " the algorithm immediately terminates with a success message."
468 " On the other hand, if the algorithm encounters 'acceptable_iter'"
469 " many iterations in a row that are considered 'acceptable', it will"
470 " terminate before the desired convergence tolerance is met. This is"
471 " useful in cases where the algorithm might not be able to achieve the"
472 "'desired' level of accuracy. The valid range for this real option is "
473 "0 < acceptable_tol < +inf and its default value is 1e-06"
474 }, 1e-6, 0, 1.e50}
475 );
476
477 slv_param_int(parameters,IPOPT_PARAM_ACCEPTABLE_ITER
478 ,(SlvParameterInitInt){{"acceptable_iter"
479 ,"Num. of 'acceptable' iters before triggering stop.",3
480 ,"If the algorithm encounters this many successive 'acceptable' "
481 "iterates (see 'acceptable_tol'), it terminates, assuming that "
482 "the problem has been solved to best possible accuracy given round-off."
483 " If it is set to zero, this heuristic is disabled. The valid range for"
484 " this integer option is 0 < acceptable_iter < +inf and its default "
485 "value is 15."
486 }, 15, 0, 100000000}
487 );
488
489
490 /** Linear Solver Options*/
491
492 /* see http://www.coin-or.org/Ipopt/documentation/node139.html */
493 slv_param_char(parameters,IPOPT_PARAM_LINEAR_SOLVER
494 ,(SlvParameterInitChar){{"linear_solver"
495 ,"Linear solver used for step computations.",4
496 ,"Determines which linear algebra package is to be used for the"
497 " solution of the augmented linear system (for obtaining the search"
498 " directions). Note, the code must have been compiled with the"
499 " linear solver you want to choose. Depending on your Ipopt"
500 " installation, not all options are available. The default value"
501 " for this string option is 'ma27'."
502 " Available options *may* include: ma27, ma57, pardiso, wsmp,"
503 " mumps, custom."
504 }, "mumps"}, (char *[]){
505 "ma27","ma57","pardiso","wsmp","mumps","custom",NULL
506 }
507 );
508
509
510 /** Barrier Parameter Options*/
511
512 slv_param_char(parameters,IPOPT_PARAM_MU_STRATEGY
513 ,(SlvParameterInitChar){{"mu_strategy"
514 ,"Update strategy for barrier parameter",5
515 ,"Determines which barrier parameter update strategy is to be used."
516 " 'monotone' is the monotone (Fiacco-McCormick) strategy;"
517 " 'adaptive' is the adaptive update strategy."
518 }, "monotone"}, (char *[]){
519 "monotone","adaptive",NULL
520 }
521 );
522
523 /** Derivative Test Options */
524
525 slv_param_char(parameters,IPOPT_PARAM_DERIVATIVE_TEST
526 ,(SlvParameterInitChar){{"derivative_test"
527 ,"Use Derivative Checker?",6
528 ,"A finite-difference derivative checker is provided by IPOPT, which"
529 " will check Jacobian and gradient functions ('first-order') or"
530 " all first-order derivatives as well as the Hessian matrix"
531 " ('second-order'). The default is to perform no checks ('none')."
532 }, "none"}, (char *[]){
533 "none","first-order","second-order",NULL
534 }
535 );
536
537 /** Quasi-Newton Options*/
538
539 slv_param_char(parameters,IPOPT_PARAM_HESS_APPROX
540 ,(SlvParameterInitChar){{"hessian_approximation"
541 ,"Hessian calculation method",7
542 ,"Use either an exact Hessian matrix based on symbolic derivatives"
543 " computed from the equations in the model ('exact'), or else use"
544 " a limited-memory quasi-Newton approximation ('limited-memory')."
545 " The default is 'limited-memory'."
546 }, "exact"}, (char *[]){
547 "exact","limited-memory",NULL
548 }
549 );
550
551
552 asc_assert(parameters->num_parms==IPOPT_PARAMS);
553
554 return 1;
555 }
556
557 static void ipopt_get_parameters(slv_system_t server, SlvClientToken asys
558 , slv_parameters_t *parameters
559 ){
560 IpoptSystem *sys;
561 UNUSED_PARAMETER(server);
562
563 sys = SYS(asys);
564 //if(check_system(sys)) return;
565 mem_copy_cast(&(sys->p),parameters,sizeof(slv_parameters_t));
566 }
567
568
569 static void ipopt_set_parameters(slv_system_t server, SlvClientToken asys
570 ,slv_parameters_t *parameters
571 ){
572 IpoptSystem *sys;
573 UNUSED_PARAMETER(server);
574
575 sys = SYS(asys);
576 //if (check_system(sys)) return;
577 mem_copy_cast(parameters,&(sys->p),sizeof(slv_parameters_t));
578 }
579
580 /*------------------------------------------------------------------------------
581 EVALUATION AND RESULT HOOK FUNCTIONS
582 */
583
584 /**
585 update the model with new 'x' vector.
586 @return 0 on success.
587 */
588 int ipopt_update_model(IpoptSystem *sys, const double *x){
589 unsigned j;
590
591 //CONSOLE_DEBUG("Updating Model ...");
592
593 asc_assert(sys);
594 asc_assert(sys->vlist);
595
596 /* FIXME do we need to update any other stuff? */
597 for(j = 0; j < sys->n; ++j){
598 //CONSOLE_DEBUG("value of var[%d] = %g", j, x[j]);
599 asc_assert(!isnan(x[j]));
600 var_set_value(sys->vlist[j], x[j]);
601 }
602
603 return 0;
604 }
605
606 /** Function to evaluate the objective function f(x).
607 @return 1 on success, 0 on failure
608
609 @param n (in), the number of variables in the problem (dimension of 'x').
610 @param x (in), the values for the primal variables, 'x' , at which 'f(x)' is to be evaluated.
611 @param new_x (in), false if any evaluation method was previously called with the same values in 'x', true otherwise.
612 @param obj_value (out) the value of the objective function ('f(x)').
613 */
614 Bool ipopt_eval_f(Index n, Number *x, Bool new_x, Number *obj_value, void *user_data){
615 IpoptSystem *sys;
616 sys = SYS(user_data);
617 int res;
618
619 //CONSOLE_DEBUG("ipopt_eval_f");
620
621 asc_assert(n==sys->n);
622 asc_assert(sys->obj!=NULL);
623
624 if(new_x){
625 res = ipopt_update_model(sys,x);
626 if(res)return 0; /* fail model update */
627 }
628
629 sys->calc_ok = TRUE;
630
631
632 /* char *relname;
633 relname = rel_make_name(sys->slv,sys->obj);
634 //CONSOLE_DEBUG("%s", relname);
635 ascfree(relname);*/
636 *obj_value = relman_eval(sys->obj,&(sys->calc_ok),SLV_PARAM_BOOL(&(sys->p),ASCEND_PARAM_SAFEEVAL));
637 asc_assert(!isnan(*obj_value));
638 //CONSOLE_DEBUG("sys->obj_value = %g",*obj_value);
639 //CONSOLE_DEBUG("done ipopt_eval_f");
640 return sys->calc_ok;
641 }
642
643 /**
644 @return 1 on success
645 */
646 Bool ipopt_eval_grad_f(Index n, Number* x, Bool new_x, Number* grad_f, void *user_data){
647 IpoptSystem *sys;
648 int j, res, len;
649 int count;
650 double *derivatives;
651 int *variables;
652
653 sys = SYS(user_data);
654
655 //CONSOLE_DEBUG("ipopt_eval_grad_f");
656
657 asc_assert(n==sys->n);
658 asc_assert(sys->obj);
659 asc_assert(sys->slv);
660
661 if(new_x){
662 res = ipopt_update_model(sys,x);
663 if(res)return 0; /* fail model update */
664 }
665
666
667 /* evaluate grad_f(x) somehow */
668 for(j=0; j<n; ++j){
669 grad_f[j] = 0;
670 }
671
672 len = rel_n_incidences(sys->obj);
673 variables = ASC_NEW_ARRAY_CLEAR(int,len);
674 derivatives = ASC_NEW_ARRAY_CLEAR(double,len);
675 /** @todo Check if memory allocation was successful and flag error if otherwise */
676 //CONSOLE_DEBUG("Length of incidences: %d",len);
677 //CONSOLE_DEBUG("allocated variables,derivatives");
678
679 /*relman_diff2(
680 sys->obj,&vfilter,derivatives,variables
681 , &count,SLV_PARAM_BOOL(&(sys->p),ASCEND_PARAM_SAFEEVAL)
682 );*/
683
684 relman_diff2_rev(
685 sys->obj,&(sys->vfilt),derivatives,variables
686 , &count,SLV_PARAM_BOOL(&(sys->p),ASCEND_PARAM_SAFEEVAL)
687 );
688
689
690 for(j=0; j<len; ++j){
691 //asc_assert(!isnan(derivatives[j]));
692 grad_f[variables[j]] = derivatives[j];
693 char *tmp = var_make_name(sys->slv, sys->vlist[variables[j]]);
694 //CONSOLE_DEBUG("var %d ('%s'): varindex = %d, x = %g, df/dx = %f", j, tmp, variables[j], var_value(sys->vlist[variables[j]]), derivatives[j]);
695 ASC_FREE(tmp);
696 }
697
698 if(variables)ASC_FREE(variables);
699 if(derivatives)ASC_FREE(derivatives);
700
701 //CONSOLE_DEBUG("done ipopt_eval_grad_f");
702 return 1; /* success, presumably */
703 }
704
705 Bool ipopt_eval_g(Index n, Number* x, Bool new_x, Index m, Number *g, void *user_data){
706 IpoptSystem *sys;
707 sys = SYS(user_data);
708 int i, res;
709 struct rel_relation *rel;
710 int calc_ok = 1;
711
712 //CONSOLE_DEBUG("ipopt_eval_g (n=%d, m=%d)",sys->n, sys->m);
713
714 asc_assert(n==sys->n);
715 asc_assert(m==sys->m);
716
717 if(new_x){
718 res = ipopt_update_model(sys,x);
719 if(res)return 0; /* fail model update */
720 }
721
722 for(i=0;i<m;++i){
723 //CONSOLE_DEBUG("rel %d: %s.",i,(sys->rlist[i] == sys->obj ? "OBJECTIVE" : "constraint")); //minor fix was rlist[0] -- MNM
724 }
725
726 /** @todo constraint rels are all relations except the objective rel. do we need to sort the objective to the end? */
727 for(i=0; i<m; ++i){
728 rel = sys->rlist[i];
729 asc_assert(rel!=NULL);
730 //if(rel == sys->obj) continue; /* I think this completes the function for the time being */
731 g[i] = relman_eval(rel, &calc_ok,SLV_PARAM_BOOL(&(sys->p),ASCEND_PARAM_SAFEEVAL));
732 asc_assert(!isnan(g[i]));
733 //CONSOLE_DEBUG("g[%d] = %f",i,g[i]);
734 }
735
736 return calc_ok; /* fail: not yet implemented */
737 }
738
739 Bool ipopt_eval_jac_g(Index n, Number* x, Bool new_x, Index m
740 , Index nele_jac, Index* iRow, Index *jCol, Number* values
741 , void *user_data
742 ){
743 IpoptSystem *sys;
744 sys = SYS(user_data);
745 int i,res,j,k,len,count;
746 struct var_variable **incidence_list;
747 int *variables;
748 double *derivatives;
749
750 //CONSOLE_DEBUG("ipopt_eval_jac_g... nnzJ = %d",sys->nnzJ);
751 //CONSOLE_DEBUG("ipopt_eval_jac_g... n = %d",sys->n);
752 //CONSOLE_DEBUG("ipopt_eval_jac_g... m = %d",sys->m);
753
754 asc_assert(sys!=NULL);
755 asc_assert(n==sys->n);
756 asc_assert(nele_jac==sys->nnzJ);
757 asc_assert(m==sys->m);
758
759 if(new_x){
760 res = ipopt_update_model(sys,x);
761 if(res)return 0; /* fail model update */
762 }
763
764 if(values == NULL){
765 CONSOLE_DEBUG("sparsity structure requested");
766 k=0;
767 for(i=0; i<m;++i){
768 /* looping through rows, one per relation */
769 if(rel_apply_filter(sys->rlist[i], &(sys->rfilt))){
770 incidence_list = (struct var_variable**) rel_incidence_list(sys->rlist[i]);
771 len=rel_n_incidences(sys->rlist[i]);
772 for(j=0;j<len;j++){
773 /* looping through incident variables in current relation */
774 if(var_apply_filter(incidence_list[j], &(sys->vfilt))){
775 CONSOLE_DEBUG("Non-zero #%d at [%d,%d]",k, i,incidence_list[j]->sindex);
776
777 /* valgrind says invalid write of size 4 here... */
778 iRow[k]=i; // should i use sindex of row here or is this ok?
779 jCol[k++]=incidence_list[j]->sindex;
780 }
781 }
782 }else{
783 CONSOLE_DEBUG("Filter removes relation %d",i);
784 }
785 }
786 CONSOLE_DEBUG("Found %d non-zero elements in jacobian", k);
787 }else{
788 //CONSOLE_DEBUG("Calculating jacobian...");
789 k=0;
790 /** @TODO make use of some temporary allocated memory for these arrays... */
791 variables = ASC_NEW_ARRAY(int,n);
792 derivatives = ASC_NEW_ARRAY_CLEAR(double,n);
793 for(i=0; i<m;++i){
794 if(rel_apply_filter(sys->rlist[i], &(sys->rfilt))){
795 incidence_list = (struct var_variable**) rel_incidence_list(sys->rlist[i]);
796 len = rel_n_incidences(sys->rlist[i]);
797
798 #if 0
799 relman_diff2(sys->rlist[i],&(sys->vfilt),derivatives,variables
800 ,&count,SLV_PARAM_BOOL(&(sys->p),ASCEND_PARAM_SAFEEVAL)
801 );
802 #else
803 relman_diff2_rev(sys->rlist[i], &(sys->vfilt), derivatives
804 ,variables, &count, SLV_PARAM_BOOL(&(sys->p),ASCEND_PARAM_SAFEEVAL)
805 );
806 #endif
807
808 for(j=0;j<count;j++){ /* loop through only the returned (filtered) incidences, not all of them */
809 asc_assert(k < sys->nnzJ);
810 //CONSOLE_DEBUG("Recording values[%d] = derivatives[%d]",k,j);
811 asc_assert(!isnan(derivatives[j]));
812 values[k++] = derivatives[j];
813 }
814 }
815 }
816 if(variables)ASC_FREE(variables);
817 if(derivatives)ASC_FREE(derivatives);
818 //CONSOLE_DEBUG("Filled in values of Jacobian");
819 }
820 //CONSOLE_DEBUG("done ipopt_eval_jac_g");
821 return TRUE;
822 }
823
824 Bool ipopt_eval_h(Index n, Number* x, Bool new_x
825 , Number obj_factor, Index m, Number* lambda
826 , Bool new_lambda, Index nele_hess, Index* iRow
827 , Index* jCol, Number* values
828 , void *user_data
829 ){
830 IpoptSystem *sys;
831 sys = SYS(user_data);
832
833 int res,count;
834
835 struct var_variable **incidence_list;
836
837 hessian_mtx *hess_matrix;
838
839 unsigned long i;
840
841 Index row;
842 Index col;
843 Index idx;
844
845 //CONSOLE_DEBUG("IN FUNCTION ipopt_eval_h");
846 //CONSOLE_DEBUG("nnzH = %d",sys->nnzH);
847 //CONSOLE_DEBUG("n = %d, m = %d",sys->n, sys->m);
848
849 asc_assert(sys!=NULL);
850 asc_assert(n==sys->n);
851 asc_assert(nele_hess==sys->nnzH);
852
853 if(new_x){
854 res = ipopt_update_model(sys,x);
855 if(res)return FALSE; /* fail model update */
856 }
857
858 if(values == NULL){
859 asc_assert(iRow !=NULL && jCol != NULL);
860
861 CONSOLE_DEBUG("Determining sparsity structure of the hessian of the lagrangian");
862
863 /* identify the sparsity structure of the Hessian (note: only the lower-
864 left part is required by IPOPT , because the Hessian is symmetric) */
865 //CONSOLE_DEBUG("Analysing of Hessian matrix sparsity structure not implemented");
866 //CONSOLE_DEBUG("Dense Hessian Calculations being performed");
867
868 idx = 0;
869
870 for (row = 0; row < n; row++) {
871 for (col = 0; col <= row; col++) {
872 iRow[idx] = row;
873 jCol[idx] = col;
874 idx++;
875 }
876 }
877 asc_assert(idx == nele_hess);
878
879 CONSOLE_DEBUG("Done with sparsity calc, there are %d elements",idx);
880 }
881 else{
882 asc_assert(jCol==NULL && iRow==NULL);
883 asc_assert(lambda!=NULL);
884
885 /** Array of LT matrix */
886 hess_matrix = Hessian_Mtx_create(Lower,n);
887
888 //CONSOLE_DEBUG("Order of Hessian MATRIX [%d x %d]",n,n);
889
890 /** Correction for objective function **/
891 //CONSOLE_DEBUG("Correction for Objective Relation underway");
892 relman_hess(sys->obj,&(sys->vfilt),hess_matrix,&count,n,SLV_PARAM_BOOL(&(sys->p),ASCEND_PARAM_SAFEEVAL));
893
894 idx = 0;
895
896 for (row = 0; row < n; row++) {
897 for (col = 0; col <= row; col++) {
898 values[idx] = Hessian_Mtx_get_element(hess_matrix,row,col) * (obj_factor);
899 idx++;
900 }
901 }
902 asc_assert(idx == nele_hess);
903
904
905 /** Correction for m-relations **/
906
907
908 for(i=0; i<m; i++){
909 /** @TODO Initialize the Hess Matrix Elements to zero */
910 Hessian_Mtx_clear(hess_matrix);
911
912 incidence_list = (struct var_variable**) rel_incidence_list(sys->rlist[i]);
913 if(incidence_list!=NULL){
914 //CONSOLE_DEBUG("Correction for Constraint Relation [%lu] underway",i);
915 relman_hess(sys->rlist[i],&(sys->vfilt),hess_matrix,&count,n,SLV_PARAM_BOOL(&(sys->p),ASCEND_PARAM_SAFEEVAL));
916
917 idx=0;
918
919 for (row = 0; row < n; row++) {
920 for (col = 0; col <= row; col++) {
921 values[idx] += Hessian_Mtx_get_element(hess_matrix,row,col) * (lambda[i]);
922 idx++;
923 }
924 }
925 asc_assert(idx == nele_hess);
926
927 }
928 else{
929 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Unused Relation???");
930 Hessian_Mtx_destroy(hess_matrix);
931 return FALSE; //I'm not sure about the action to take.
932 }
933 }
934
935 //CONSOLE_DEBUG("Hessian Matrix evaluation successful");
936
937 Hessian_Mtx_destroy(hess_matrix);
938
939 /* evaluate the Hessian matrix */
940 //CONSOLE_DEBUG("Evaluation of Hessian matrix Completed");
941 }
942
943 return TRUE; /* fail: not yet implemented */
944 }
945
946 /*------------------------------------------------------------------------------
947 SOLVE ROUTINES
948 */
949
950 static int ipopt_presolve(slv_system_t server, SlvClientToken asys){
951 IpoptSystem *sys;
952 int max, i;
953 struct var_variable *var;
954
955 //CONSOLE_DEBUG("PRESOLVE");
956
957 sys = SYS(asys);
958 ipopt_iteration_begins(sys);
959 //check_system(sys);
960
961 asc_assert(sys->vlist && sys->rlist);
962
963 /** @todo work out if matrix creation is not again needed */
964
965 slv_sort_rels_and_vars(server,&(sys->m),&(sys->n));
966 #if 0
967 /* ignore any errors here; if it fails, we may just have a single objective function and no constraining relations */
968 if(-1 == sys->n){
969 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to find any optimisable vars");
970 return -4;
971 }
972 if(-1 == sys->m){
973 sys->m = 0; /* no relations found, but that's OK if there's an objective? */
974 }
975 if(-1 == sys->m)sys->m = 0;
976 if(-1 == sys->n)sys->n = 0;
977 #endif
978
979 CONSOLE_DEBUG("Got %d rels and %d vars",sys->m, sys->n);
980
981 #if 1
982 /* count the number of optimisation variables */
983 sys->n = 0;
984 for(i = 0; i < sys->vtot; i++){
985 var = sys->vlist[i];
986 if(var_apply_filter(var,&(sys->vfilt))){
987 sys->n++;
988 }
989 }
990 #endif
991
992 /* set all relations as being 'unsatisfied' to start with... */
993 for(i=0; i < sys->rtot; ++i){
994 rel_set_satisfied(sys->rlist[i],FALSE);
995 }
996
997 sys->obj = slv_get_obj_relation(server); /*may have changed objective*/
998
999
1000 if(!sys->obj){
1001 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No objective function was specified");
1002 return -3;
1003 }
1004 //CONSOLE_DEBUG("got objective rel %p",sys->obj);
1005 /* @todo check if old_obj == obj ? */
1006
1007 #if 1
1008 /* TODO are there cases where these should be different: answer: NO. they are always the same -- JP */
1009 sys->m = sys->rtot;
1010 #endif
1011
1012 //CONSOLE_DEBUG("Numbers of constraints = %d",sys->m);
1013
1014 /** @todo we need to move the objective relation to the end of the list */
1015
1016 /*for(i=0;i<sys->rtot-1;++i){
1017 //CONSOLE_DEBUG("%d",i);
1018 if(sys->rlist[i] == sys->obj)
1019 //CONSOLE_DEBUG("<-------------------------------This Check Works------------------------>");
1020
1021 }*/
1022
1023 //CONSOLE_DEBUG("got objective rel %p",sys->obj);
1024
1025 /* calculate nnz for hessian matrix @todo FIXME */
1026
1027 if(strcmp(SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_HESS_APPROX),"exact")==0){
1028 /** @todo fix rtot to be 'm' instead */
1029 sys->nnzH = ((sys->n)*((sys->n)+1))/2; //dense Hessian count
1030 }else{
1031 //CONSOLE_DEBUG("Skipping relman_hessian_count as hessian method is not exact.");
1032 //sys->nnzH = sys->n * sys->m;
1033 }
1034
1035 /* need to provide sparsity structure for hessian matrix? */
1036
1037 #if 0
1038 /** @SEE http://www.coin-or.org/Ipopt/documentation/node37.html */
1039 ipopt_eval_h(number_of_variables, NULL/*x at which to evaluate*/, TRUE /* new x */
1040 , 1.0/*obj_factor*/, number_of_constraints, lambda/* values of the constraint multipliers */
1041 , TRUE /* new lambda */, 0 /* number of nonzero elements in the Hessian */, Index* iRow
1042 , Index* jCol, Number* values
1043 , void *user_data
1044 );
1045 #endif
1046
1047
1048 max = relman_obj_direction(sys->obj);
1049 if(max==-1){
1050 //CONSOLE_DEBUG("this is a MINIMIZE problem");
1051 }else{
1052 //CONSOLE_DEBUG("this is a MAXIMIZE problem");
1053 }
1054
1055 //CONSOLE_DEBUG("got %d constraints and %d vars in system", sys->m, sys->n);
1056 /* calculate number of non-zeros in the Jacobian matrix for the constraint equations */
1057
1058 /* @todo make sure objective rel moved to end */
1059
1060 CONSOLE_DEBUG("About to call relman_jacobian_count");
1061 sys->nnzJ = relman_jacobian_count(sys->rlist, sys->m, &(sys->vfilt), &(sys->rfilt), &max);
1062 /*sys->nnzJ=0;
1063 for(i=0;i<sys->m;++i){
1064 sys->nnzJ += rel_n_incidences(sys->rlist[i]);
1065 }*/
1066
1067 CONSOLE_DEBUG("got %d non-zeros in constraint Jacobian", sys->nnzJ);
1068
1069 /* need to provide sparsity structure for jacobian? */
1070
1071
1072
1073 #if 0
1074 if(sys->presolved > 0) { /* system has been presolved before */
1075 if(!conopt_dof_changed(sys) /*no changes in fixed or included flags*/
1076 && sys->p.partition == sys->J.old_partition
1077 && sys->obj == sys->old_obj
1078 ){
1079 matrix_creation_needed = 0;
1080 CONOPT_//CONSOLE_DEBUG("YOU JUST AVOIDED MATRIX DESTRUCTION/CREATION");
1081 }
1082 }
1083 #endif
1084
1085 #if 0
1086 // check all this...
1087
1088 sys->presolved = 1; /* full presolve recognized here */
1089 sys->resolve = 0; /* initialize resolve flag */
1090
1091 sys->J.old_partition = sys->p.partition;
1092 sys->old_obj = sys->obj;
1093
1094 slv_sort_rels_and_vars(server,&(sys->con.m),&(sys->con.n));
1095 CONOPT_//CONSOLE_DEBUG("FOUND %d CONSTRAINTS AND %d VARS",sys->con.m,sys->con.n);
1096 if (sys->obj != NULL) {
1097 CONOPT_//CONSOLE_DEBUG("ADDING OBJECT AS A ROW");
1098 sys->con.m++; /* treat objective as a row */
1099 }
1100
1101 cntvect = ASC_NEW_ARRAY(int,COIDEF_Size());
1102 COIDEF_Ini(cntvect);
1103 sys->con.cntvect = cntvect;
1104 CONOPT_//CONSOLE_DEBUG("NUMBER OF CONSTRAINTS = %d",sys->con.m);
1105 COIDEF_NumVar(cntvect, &(sys->con.n));
1106 COIDEF_NumCon(cntvect, &(sys->con.m));
1107 sys->con.nz = num_jacobian_nonzeros(sys, &(sys->con.maxrow));
1108 COIDEF_NumNZ(cntvect, &(sys->con.nz));
1109 COIDEF_NumNlNz(cntvect, &(sys->con.nz));
1110
1111 sys->con.base = 0;
1112 COIDEF_Base(cntvect,&(sys->con.base));
1113 COIDEF_ErrLim(cntvect, &(DOMLIM));
1114 COIDEF_ItLim(cntvect, &(ITER_LIMIT));
1115
1116 if(sys->obj!=NULL){
1117 sys->con.optdir = relman_obj_direction(sys->obj);
1118 sys->con.objcon = sys->con.m - 1; /* objective will be last row */
1119 CONOPT_//CONSOLE_DEBUG("SETTING OBJECTIVE CONSTRAINT TO BE %d",sys->con.objcon);
1120 }else{
1121 sys->con.optdir = 0;
1122 sys->con.objcon = 0;
1123 }
1124 COIDEF_OptDir(cntvect, &(sys->con.optdir));
1125 COIDEF_ObjCon(cntvect, &(sys->con.objcon));
1126
1127 temp = 0;
1128 COIDEF_StdOut(cntvect, &temp);
1129
1130 int debugfv = 1;
1131 COIDEF_DebugFV(cntvect, &debugfv);
1132
1133 destroy_vectors(sys);
1134 destroy_matrices(sys);
1135 create_matrices(server,sys);
1136 create_vectors(sys);
1137
1138 sys->s.block.current_reordered_block = -2;
1139 }
1140
1141 //...
1142
1143 if( matrix_creation_needed ) {
1144 destroy_array(sys->s.cost);
1145 sys->s.cost = create_zero_array(sys->s.costsize,struct slv_block_cost);
1146 for( ind = 0; ind < sys->s.costsize; ++ind ) {
1147 sys->s.cost[ind].reorder_method = -1;
1148 }
1149 } else {
1150 reset_cost(sys->s.cost,sys->s.costsize);
1151 }
1152
1153 #endif
1154
1155 /* Reset status */
1156 sys->s.iteration = 0;
1157 sys->s.cpu_elapsed = 0.0;
1158 sys->s.converged = sys->s.diverged = sys->s.inconsistent = FALSE;
1159 sys->s.block.previous_total_size = 0;
1160 sys->s.costsize = 1+sys->s.block.number_of;
1161
1162
1163 /* set to go to first unconverged block */
1164 sys->s.block.current_block = -1;
1165 sys->s.block.current_size = 0;
1166 sys->s.calc_ok = TRUE;
1167 sys->s.block.iteration = 0;
1168 sys->obj_val = MAXDOUBLE/2000.0;
1169 //CONSOLE_DEBUG("sys->obj_val=%g",sys->obj_val);
1170 update_status(sys);
1171
1172 ipopt_iteration_ends(sys);
1173
1174 //CONSOLE_DEBUG("Reset status");
1175
1176 /* sys->s.cost[sys->s.block.number_of].time=sys->s.cpu_elapsed; */
1177
1178 //ERROR_REPORTER_HERE(ASC_USER_SUCCESS,"presolve completed");
1179 return 0;
1180 }
1181
1182
1183 static int ipopt_solve(slv_system_t server, SlvClientToken asys){
1184 IpoptSystem *sys;
1185 UNUSED_PARAMETER(server);
1186 enum ApplicationReturnStatus status;
1187 int ret, i, j;
1188 struct var_variable *var;
1189 enum rel_enum type_of_rel;
1190 sys = SYS(asys);
1191
1192 double *x, *x_L, *x_U, *g_L, *g_U, *mult_x_L, *mult_x_U;
1193
1194 CONSOLE_DEBUG("SOLVING: sys->n = %d, sys->m = %d...",sys->n,sys->m);
1195 asc_assert(sys->n!=-1);
1196
1197 /* set the number of variables and allocate space for the bounds */
1198 x_L = ASC_NEW_ARRAY(Number,sys->n);
1199 x_U = ASC_NEW_ARRAY(Number,sys->n);
1200
1201 //CONSOLE_DEBUG("SETTING BOUNDS...");
1202
1203 /* @todo set the values for the variable bounds */
1204 int jj = 0;
1205 for(j = 0; j < sys->vtot; j++){
1206 //CONSOLE_DEBUG("j = %d, vtot = %d, vlist = %p",j,sys->vtot,sys->vlist);
1207 var = sys->vlist[j];
1208 if(var_apply_filter(var,&(sys->vfilt))){
1209 //CONSOLE_DEBUG("setting x_L[%d] = %e",jj,var_lower_bound(var));
1210 assert(jj<sys->n);
1211 x_L[jj] = var_lower_bound(var);
1212 //CONSOLE_DEBUG("setting x_U[%d] = %e",jj,var_upper_bound(var));
1213 x_U[jj] = var_upper_bound(var);
1214 jj++;
1215 }
1216 }
1217
1218 //CONSOLE_DEBUG("jj = %d, sys->n = %d", jj, sys->n);
1219 assert(jj==sys->n);
1220
1221 /** @todo set bounds on the constraints? */
1222 /* is it possible to identify f(x)<a; f(x) >b and fold them into one? */
1223 /* then find the constant parts and make then g_L or g_U accordingly */
1224 /* what to do about other bounds? */
1225 /* set the number of variables and allocate space for the bounds */
1226 g_L = ASC_NEW_ARRAY(Number,sys->m);
1227 g_U = ASC_NEW_ARRAY(Number,sys->m);
1228 //CONSOLE_DEBUG("Allocated arrays for bounds of relations");
1229 if(g_L!=NULL && g_U!=NULL)
1230 for(j = 0; j < sys->m; j++){
1231 type_of_rel = rel_relop(sys->rlist[j]);
1232 if (type_of_rel == e_rel_less || type_of_rel == e_rel_lesseq){
1233 g_L[j] = -2.0e19; //refer to IPOPT Manual "The C Interface"
1234 g_U[j] = 0;
1235 }
1236 else if (type_of_rel == e_rel_greatereq || type_of_rel == e_rel_greater){
1237 g_L[j] = 0;
1238 g_U[j] = 2.0e19; //refer to IPOPT Manual "the C Interface"
1239 }
1240 else{
1241 g_L[j] = 0;
1242 g_U[j] = 0;
1243 }
1244 //CONSOLE_DEBUG("set g_L[%d] = %e",j,g_L[j]);
1245 //CONSOLE_DEBUG("set g_U[%d] = %e",j,g_U[j]);
1246 }
1247 else
1248 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to allocate arrays for bounds of relations");
1249
1250
1251 //CONSOLE_DEBUG("CREATING PROBLEM...");
1252
1253 /* create the IpoptProblem */
1254 //CONSOLE_DEBUG("n = %d, m = %d, nnzJ = %d, nnzH = %d",sys->n, sys->m, sys->nnzJ, sys->nnzH);
1255 sys->nlp = CreateIpoptProblem(sys->n, x_L, x_U, sys->m, g_L, g_U, sys->nnzJ, sys->nnzH, 0/*index style=C*/,
1256 &ipopt_eval_f, &ipopt_eval_g, &ipopt_eval_grad_f,
1257 &ipopt_eval_jac_g, &ipopt_eval_h
1258 );
1259
1260 //CONSOLE_DEBUG("FREEING INTERNAL STUFF");
1261
1262 /* We can free the memory now - the values for the bounds have been
1263 copied internally in CreateIpoptProblem */
1264 ASC_FREE(x_L);
1265 ASC_FREE(x_U);
1266 ASC_FREE(g_L);
1267 ASC_FREE(g_U);
1268
1269 //CONSOLE_DEBUG("SETTING OPTIONS...");
1270 /* set some options */
1271 /** OUTPUT OPTIONS */
1272 AddIpoptIntOption(sys->nlp, "print_level", SLV_PARAM_INT(&(sys->p),IPOPT_PARAM_PRINT_LEVEL));
1273 AddIpoptStrOption(sys->nlp, "print_user_options", SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_PRINT_USER_OPTIONS));
1274 /** TERMINATION OPTIONS */
1275 AddIpoptIntOption(sys->nlp, "max_iter", SLV_PARAM_INT(&(sys->p),IPOPT_PARAM_MAX_ITER));
1276 AddIpoptNumOption(sys->nlp, "tol", SLV_PARAM_REAL(&(sys->p),IPOPT_PARAM_TOL));
1277 AddIpoptNumOption(sys->nlp, "max_cpu_time", SLV_PARAM_REAL(&(sys->p),IPOPT_PARAM_MAX_CPU_TIME));
1278 AddIpoptNumOption(sys->nlp, "diverging_iterates_tol", SLV_PARAM_REAL(&(sys->p),IPOPT_PARAM_DIVERGING_ITERATES_TOL));
1279 AddIpoptNumOption(sys->nlp, "dual_inf_tol", SLV_PARAM_REAL(&(sys->p),IPOPT_PARAM_DUAL_INFEASIBILITY_TOL));
1280 AddIpoptNumOption(sys->nlp, "constr_viol_tol", SLV_PARAM_REAL(&(sys->p),IPOPT_PARAM_CONSTR_VIOL_TOL));
1281 AddIpoptNumOption(sys->nlp, "acceptable_tol", SLV_PARAM_REAL(&(sys->p),IPOPT_PARAM_ACCEPTABLE_TOL));
1282 AddIpoptIntOption(sys->nlp, "acceptable_iter", SLV_PARAM_INT(&(sys->p),IPOPT_PARAM_ACCEPTABLE_ITER));
1283 /** BARRIER PARAMETER OPTIONS */
1284 AddIpoptStrOption(sys->nlp, "mu_strategy", SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_MU_STRATEGY));
1285 /** DERIVATIVE TEST OPTIONS */
1286 AddIpoptStrOption(sys->nlp, "derivative_test", SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_DERIVATIVE_TEST));
1287 /** QUASI-NEWTON OPTIONS */
1288 AddIpoptStrOption(sys->nlp, "hessian_approximation", SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_HESS_APPROX));
1289 /** LINEAR SOLVER OPTIONS */
1290 AddIpoptStrOption(sys->nlp, "linear_solver", SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_LINEAR_SOLVER));
1291
1292
1293 //CONSOLE_DEBUG("Hessian method: %s",SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_HESS_APPROX));
1294
1295 //CONSOLE_DEBUG("number of vars n = %d, number of rels m = %d",sys->n, sys->m);
1296
1297 /* initial values */
1298 x = ASC_NEW_ARRAY(Number, sys->n);
1299 /*setting initial values here.*/
1300 //CONSOLE_DEBUG("Setting starting values for free variables.");
1301 for(i=0;i<sys->n;++i){
1302 //CONSOLE_DEBUG("set x[%d] = %g",i,var_value(sys->vlist[i])); // need to set the default values
1303 x[i]=var_value(sys->vlist[i]);
1304 }
1305 /** @todo get values of 'x' from the model */
1306
1307 /* allocate space to store the bound multipliers at the solution */
1308 mult_x_L = ASC_NEW_ARRAY(Number, sys->n);
1309 mult_x_U = ASC_NEW_ARRAY(Number, sys->n);
1310
1311 //CONSOLE_DEBUG("Calling IpoptSolve...");
1312
1313 //CONSOLE_DEBUG("sys->objval = %g", sys->obj_val);
1314
1315 /* solve the problem */
1316 status = IpoptSolve(sys->nlp, x, NULL, &sys->obj_val, NULL, mult_x_L, mult_x_U, (void*)sys);
1317
1318 //CONSOLE_DEBUG("Done IpoptSolve...");
1319
1320 /** @todo update the sys->s.xxxxx flags based on value of 'status' */
1321
1322 ret = 1; /* default case is failure */
1323 switch(status){
1324 case Solve_Succeeded:
1325 sys->s.converged = TRUE;
1326
1327 sys->s.block.current_block = -1; //is this 1??
1328 sys->s.cost = ASC_NEW_ARRAY(struct slv_block_cost,1);
1329 sys->s.cost->size=sys->s.block.current_size=sys->n;
1330 sys->s.cost->iterations=sys->s.block.iteration;
1331 sys->s.cost->funcs=sys->s.block.funcs;
1332 sys->s.cost->jacs=sys->s.block.jacs;
1333 sys->s.cost->time=sys->s.block.cpu_elapsed;
1334 sys->s.cost->functime=sys->s.block.functime;
1335 sys->s.cost->jactime=sys->s.block.jactime;
1336
1337
1338 //CONSOLE_DEBUG("Solution of the primal variables, x");
1339 for (i=0; i<sys->n; i++) {
1340 //CONSOLE_DEBUG(" x[%d] = %e\n", i, x[i]);
1341 }
1342
1343 //CONSOLE_DEBUG("Solution of the bound multipliers, z_L and z_U");
1344 for (i=0; i<sys->n; i++) {
1345 //CONSOLE_DEBUG(" z_L[%d] = %e", i, mult_x_L[i]);
1346 }
1347 for (i=0; i<sys->n; i++) {
1348 //CONSOLE_DEBUG(" z_U[%d] = %e", i, mult_x_U[i]);
1349 }
1350
1351 //CONSOLE_DEBUG("Objective value");
1352 //CONSOLE_DEBUG(" f(x*) = %e", sys->obj_val);
1353
1354 ret = 0; /* success */
1355 ipopt_iteration_ends(sys);
1356 update_status(sys);
1357
1358 break;
1359 case Search_Direction_Becomes_Too_Small:
1360 ERROR_REPORTER_HERE(ASC_USER_NOTE,"Solve direction becomes too small");
1361 break;
1362 case Feasible_Point_Found:
1363 ERROR_REPORTER_HERE(ASC_USER_NOTE,"Feasible point not found");
1364 break;
1365 case NonIpopt_Exception_Thrown:
1366 ERROR_REPORTER_HERE(ASC_USER_NOTE,"Non-IPOPT exception thrown");
1367 break;
1368 case Solved_To_Acceptable_Level:
1369 /** @todo What should be done here? */
1370 ERROR_REPORTER_HERE(ASC_USER_NOTE,"Solved to acceptable level");
1371 break;
1372 case Infeasible_Problem_Detected:
1373 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Infeasible Problem Detected");
1374 break;
1375 case Diverging_Iterates:
1376 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Diverging iterations found.");
1377 break;
1378 case User_Requested_Stop:
1379 ERROR_REPORTER_HERE(ASC_USER_WARNING,"User Requested Stop.");
1380 break;
1381 case Maximum_Iterations_Exceeded:
1382 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Maximum Iterations Exceeded.");
1383 break;
1384 case Restoration_Failed:
1385 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Restoration Failed.");
1386 break;
1387 case Error_In_Step_Computation:
1388 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Error in Step Computation.");
1389 break;
1390 case Maximum_CpuTime_Exceeded:
1391 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Maximum CPU Time exceeded.");
1392 break;
1393 case Not_Enough_Degrees_Of_Freedom:
1394 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Not enough degrees of freedom.");
1395 break;
1396 case Invalid_Problem_Definition:
1397 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Invalid problem definition.");
1398 break;
1399 case Invalid_Option:
1400 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Invalid Option.");
1401 break;
1402 case Invalid_Number_Detected:
1403 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Invalid Number Detected.");
1404 break;
1405 case Unrecoverable_Exception:
1406 ERROR_REPORTER_HERE(ASC_PROG_FATAL,"Unrecoverable_Exception.");
1407 break;
1408 case Insufficient_Memory:
1409 ERROR_REPORTER_HERE(ASC_PROG_FATAL,"Insufficient Memory.");
1410 break;
1411 case Internal_Error:
1412 ERROR_REPORTER_HERE(ASC_PROG_FATAL,"Internal Error.");
1413 break;
1414 default:
1415 ERROR_REPORTER_HERE(ASC_PROG_ERROR,"Unhanded return state %d from IPOPT",status);
1416 }
1417
1418 /* free allocated memory */
1419 FreeIpoptProblem(sys->nlp);
1420 ASC_FREE(x);
1421 ASC_FREE(mult_x_L);
1422 ASC_FREE(mult_x_U);
1423
1424 return ret;
1425 }
1426
1427 /**
1428 Prepare sys for entering an iteration, increasing the iteration counts
1429 and starting the clock.
1430 */
1431 static void ipopt_iteration_begins(IpoptSystem *sys){
1432 sys->clock = tm_cpu_time();
1433 ++(sys->s.block.iteration);
1434 ++(sys->s.iteration);
1435 }
1436
1437
1438 /*
1439 Prepare sys for exiting an iteration, stopping the clock and recording
1440 the cpu time.
1441 */
1442 static void ipopt_iteration_ends(IpoptSystem *sys){
1443 double cpu_elapsed; /* elapsed this iteration */
1444
1445 cpu_elapsed = (double)(tm_cpu_time() - sys->clock);
1446 sys->s.block.cpu_elapsed += cpu_elapsed;
1447 sys->s.cpu_elapsed += cpu_elapsed;
1448 }
1449
1450
1451
1452 static int ipopt_iterate(slv_system_t server, SlvClientToken asys){
1453 //CONSOLE_DEBUG("ipopt_iterate about to call ipopt_solve...");
1454 return ipopt_solve(server,asys);
1455 }
1456
1457 static int ipopt_resolve(slv_system_t server, SlvClientToken asys){
1458 IpoptSystem *sys;
1459 sys = SYS(asys);
1460
1461 /** @todo if implementing this, use the 'warm start' thing in IPOPT */
1462
1463 /** @todo provide initial values of the 'multipliers' */
1464
1465 sys->resolve = 1; /* resolved recognized here */
1466
1467 /* Reset status */
1468 sys->s.iteration = 0;
1469 sys->s.cpu_elapsed = 0.0;
1470 sys->s.converged = sys->s.diverged = sys->s.inconsistent = FALSE;
1471 sys->s.block.previous_total_size = 0;
1472
1473 /* go to first unconverged block */
1474 sys->s.block.current_block = -1;
1475 sys->s.block.current_size = 0;
1476 sys->s.calc_ok = TRUE;
1477 sys->s.block.iteration = 0;
1478 sys->obj_val = MAXDOUBLE/2000.0;
1479
1480 update_status(sys);
1481 return 1;
1482 }
1483
1484 static const SlvFunctionsT ipopt_internals = {
1485 67
1486 ,"IPOPT"
1487 ,ipopt_create
1488 ,ipopt_destroy
1489 ,ipopt_eligible_solver
1490 ,ipopt_get_default_parameters
1491 ,ipopt_get_parameters
1492 ,ipopt_set_parameters
1493 ,ipopt_get_status
1494 ,ipopt_solve
1495 ,ipopt_presolve
1496 ,ipopt_iterate
1497 ,ipopt_resolve
1498 ,NULL
1499 ,NULL
1500 ,NULL
1501 };
1502
1503 int ipopt_register(void){
1504 return solver_register(&ipopt_internals);
1505 }

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22