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

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