/[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 2039 - (show annotations) (download) (as text)
Wed May 20 01:34:20 2009 UTC (13 years, 8 months ago) by jpye
File MIME type: text/x-csrc
File size: 29594 byte(s)
Updated formatting/file headers for set.h, mem.h.
Added #includes to top of set.h, mem.h to ensure ascConfig is not missed.
Changed default IPOPT linear solver to MUMPS.
Added simple C++ test file for CONOPT.
1 /* ASCEND modelling environment
2 Copyright (C) 2007-2008 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.
27 */
28
29 #include <ascend/utilities/config.h>
30
31 #ifndef ASC_WITH_IPOPT
32 # error "ASC_WITH_IPOPT must be defined in order to build this."
33 #endif
34
35 #include <ascend/solver/solver.h>
36
37 #include <ascend/system/calc.h>
38 #include <ascend/system/relman.h>
39 #include <ascend/system/slv_stdcalls.h>
40 #include <ascend/system/block.h>
41
42 #include <ascend/utilities/ascConfig.h>
43 #include <ascend/utilities/ascPanic.h>
44 #include <ascend/utilities/ascMalloc.h>
45 #include <ascend/utilities/ascDynaLoad.h>
46 #include <ascend/utilities/mem.h>
47 #include <ascend/utilities/ascEnvVar.h>
48 #include <ascend/general/tm_time.h>
49 #include <ascend/general/env.h>
50
51 #include <coin/IpStdCInterface.h>
52
53 ASC_DLLSPEC SolverRegisterFn ipopt_register;
54
55 /*------------------------------------------------------------------------------
56 DATA STRUCTURES AND FORWARD DEFS
57 */
58
59 /**
60 Documentation of solver options for IPOPT is at
61 http://www.coin-or.org/Ipopt/documentation/node1.html
62 */
63 enum{
64 IPOPT_PARAM_TOL
65 ,IPOPT_PARAM_LINEAR_SOLVER
66 ,IPOPT_PARAM_MAX_ITER
67 ,IPOPT_PARAM_SAFEEVAL
68 ,IPOPT_PARAM_MU_STRATEGY
69 ,IPOPT_PARAM_DERIVATIVE_TEST
70 ,IPOPT_PARAM_HESS_APPROX
71 ,IPOPT_PARAMS
72 };
73
74 #define SYS(s) ((IpoptSystem *)(s))
75
76 struct IpoptSystemStruct{
77
78 /*
79 Problem definition
80 */
81 slv_system_t slv; /* slv_system_t back-link */
82 struct rel_relation *obj; /* Objective function: NULL = none */
83 struct rel_relation *old_obj;/* Objective function: NULL = none */
84 struct var_variable **vlist; /* Variable list (NULL terminated) */
85 struct rel_relation **rlist; /* Relation list (NULL terminated) */
86
87 var_filter_t vfilt;
88 rel_filter_t rfilt;
89
90 /*
91 Solver information
92 */
93 int32 presolved; /* ? Has the system been presolved */
94 int32 resolve; /* ? Has the system been resolved */
95 slv_parameters_t p; /* Parameters */
96 slv_status_t s; /* Status (as of iteration end) */
97
98 int32 cap; /* Order of matrix/vectors */
99 int32 rank; /* Symbolic rank of problem */
100 int32 vused; /* Free and incident variables */
101 int32 vtot; /* length of varlist */
102 int32 rused; /* Included relations */
103 int32 rtot; /* length of rellist */
104 double clock; /* CPU time */
105
106 int32 calc_ok;
107 double obj_val;
108
109 void *parm_array[IPOPT_PARAMS];
110 struct slv_parameter pa[IPOPT_PARAMS];
111
112 /*
113 IPOPT DATA
114 */
115 Index n; /* number of variables */
116 Index m; /* number of constraints (excl the 'objective relation')*/
117
118 Index nnzJ; /* number of non zeros in the jacobian of the constraints */
119 Index nnzH; /* number of non-zeros in the hessian of the objective */
120
121 Number* x_L; /* lower bounds on x */
122 Number* x_U; /* upper bounds on x */
123 Number* g_L; /* lower bounds on g */
124 Number* g_U; /* upper bounds on g */
125
126 IpoptProblem nlp; /* IpoptProblem */
127
128 enum ApplicationReturnStatus status; /* Solve return code */
129 Number* x; /* starting point and solution vector */
130 Number* mult_x_L; /* lower bound multipliers at the solution */
131 Number* mult_x_U; /* upper bound multipliers at the solution */
132 Index i; /* generic counter */
133 };
134
135 typedef struct IpoptSystemStruct IpoptSystem;
136
137 static int ipopt_get_default_parameters(slv_system_t server, SlvClientToken asys, slv_parameters_t *parameters);
138
139 static void ipopt_iteration_begins(IpoptSystem *sys);
140 static void ipopt_iteration_ends(IpoptSystem *sys);
141
142 /*------------------------------------------------------------------------------
143 SYSTEM SETUP/DESTROY, STATUS AND SOLVER ELIGIBILITY
144 */
145
146 static SlvClientToken ipopt_create(slv_system_t server, int32*statusindex){
147 IpoptSystem *sys;
148
149 sys = ASC_NEW_CLEAR(IpoptSystem);
150 if(sys==NULL){
151 *statusindex = 1;
152 return sys;
153 }
154
155 sys->p.parms = sys->pa;
156 sys->p.dynamic_parms = 0;
157 ipopt_get_default_parameters(server,(SlvClientToken)sys,&(sys->p));
158 sys->p.whose = (*statusindex);
159
160 sys->presolved = 0;
161 sys->resolve = 0;
162
163 sys->n = -1;
164 sys->m = -1;
165
166 sys->s.ok = TRUE;
167 sys->s.calc_ok = TRUE;
168 sys->s.costsize = 0;
169 sys->s.cost = NULL; /*redundant, but sanity preserving */
170 sys->s.block.number_of = 1;
171 sys->s.block.current_block = 0;
172 sys->s.block.current_reordered_block = 0;
173 sys->s.block.current_size = 0;
174 sys->s.block.previous_total_size = 0;
175 sys->s.block.iteration = 0;
176 sys->s.block.funcs = 0;
177 sys->s.block.jacs = 0;
178 sys->s.block.cpu_elapsed = 0;
179 sys->s.block.functime = 0;
180 sys->s.block.jactime = 0;
181 sys->s.block.residual = 0;
182
183 sys->rfilt.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
184 sys->rfilt.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
185 sys->vfilt.matchbits = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR | VAR_FIXED);
186 sys->vfilt.matchvalue = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR);
187
188 sys->vlist = slv_get_solvers_var_list(server);
189 sys->rlist = slv_get_solvers_rel_list(server);
190
191 sys->rtot = slv_get_num_solvers_rels(server);
192 sys->vtot = slv_get_num_solvers_vars(server);
193
194 sys->obj = slv_get_obj_relation(server);
195 if(sys->vlist == NULL) {
196 ASC_FREE(sys);
197 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IPOPT called with no variables.");
198 *statusindex = -2;
199 return NULL;
200 }
201 if(sys->rlist == NULL && sys->obj == NULL) {
202 ASC_FREE(sys);
203 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IPOPT called with no relations or objective.");
204 *statusindex = -1;
205 return NULL;
206 }
207
208 /* do nothing with the objective list, pars, bounds, extrels, etc etc */
209
210 slv_check_var_initialization(server);
211 *statusindex = 0;
212 return((SlvClientToken)sys);
213 }
214
215 static int32 ipopt_destroy(slv_system_t server, SlvClientToken asys){
216 UNUSED_PARAMETER(server);
217 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ipopt_destroy not implemented");
218 return 1;
219 }
220
221
222 static int ipopt_get_status(slv_system_t server, SlvClientToken asys
223 ,slv_status_t *status
224 ){
225 IpoptSystem *sys;
226 (void)server; /* stop gcc whine about unused parameter */
227
228 sys = SYS(asys);
229 //if (check_system(sys)) return 1;
230 mem_copy_cast(&(sys->s),status,sizeof(slv_status_t));
231 return 0;
232 }
233
234 /**
235 Update the solver status. FIXME can't we get rid of this silly function
236 somehot?
237 */
238 static void update_status(IpoptSystem *sys){
239 boolean unsuccessful;
240
241 sys->s.time_limit_exceeded = FALSE; /* can't do this one with IPOPT */
242 sys->s.iteration_limit_exceeded = FALSE; /* IPOPT handles this one internally */
243
244 unsuccessful = sys->s.diverged || sys->s.inconsistent ||
245 sys->s.iteration_limit_exceeded || sys->s.time_limit_exceeded;
246
247 sys->s.ready_to_solve = !unsuccessful && !sys->s.converged;
248 sys->s.ok = !unsuccessful && sys->s.calc_ok && !sys->s.struct_singular;
249 }
250
251 static int32 ipopt_eligible_solver(slv_system_t server){
252 UNUSED_PARAMETER(server);
253
254 /// @todo check that there is a MAXIMIZE or MINIMIZE statement
255 /// @todo check that there are no discrete-valued variables
256 /// @todo check that there are no WHENs or CONDITIONALs
257 /// @todo check anything else?
258
259 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"ipopt_eligible_solver not implemented");
260 return 1;
261 }
262
263 /*------------------------------------------------------------------------------
264 SOLVER PARAMETERS
265 */
266
267 static
268 int32 ipopt_get_default_parameters(slv_system_t server, SlvClientToken asys
269 ,slv_parameters_t *parameters
270 ){
271 IpoptSystem *sys = NULL;
272 struct slv_parameter *new_parms = NULL;
273 int32 make_macros = 0;
274
275 if(server != NULL && asys != NULL) {
276 sys = SYS(asys);
277 make_macros = 1;
278 }
279
280 if(parameters->parms == NULL) {
281 new_parms = ASC_NEW_ARRAY_OR_NULL(struct slv_parameter,IPOPT_PARAMS);
282 if(new_parms == NULL) {
283 return -1;
284 }
285 parameters->parms = new_parms;
286 parameters->dynamic_parms = 1;
287 }
288
289 parameters->num_parms = 0;
290
291 slv_param_int(parameters,IPOPT_PARAM_MAX_ITER
292 ,(SlvParameterInitInt){{"max_iter"
293 ,"Maximum number of iterations",2
294 ,"The algorithm terminates with an error message if the number of iterations exceeded this number."
295 }, 3000, 0, 100000000}
296 );
297
298 /* see http://www.coin-or.org/Ipopt/documentation/node139.html */
299 slv_param_char(parameters,IPOPT_PARAM_LINEAR_SOLVER
300 ,(SlvParameterInitChar){{"linear_solver"
301 ,"Linear solver used for step computations.",2
302 ,"Determines which linear algebra package is to be used for the"
303 " solution of the augmented linear system (for obtaining the search"
304 " directions). Note, the code must have been compiled with the"
305 " linear solver you want to choose. Depending on your Ipopt"
306 " installation, not all options are available. The default value"
307 " for this string option is 'ma27'."
308 " Available options *may* include: ma27, ma57, pardiso, wsmp,"
309 " mumps, custom."
310 }, "mumps"}, (char *[]){
311 "ma27","ma57","pardiso","wsmp","mumps","custom",NULL
312 }
313 );
314
315 slv_param_real(parameters,IPOPT_PARAM_TOL
316 ,(SlvParameterInitReal){{"tol"
317 ,"Desired convergence tolerance (relative)",2
318 ,"Determines the convergence tolerance for the algorithm. The algorithm"
319 " terminates successfully, if the (scaled) NLP error becomes smaller"
320 " than this value, and if the (absolute) criteria according to "
321 " 'dual_inf_tol', 'primal_inf_tol', and 'cmpl_inf_tol' are met. (This"
322 " is epsilon_tol in Eqn. (6) in implementation paper). See also "
323 " 'acceptable_tol' as a second termination criterion. Note, some other"
324 " algorithmic features also use this quantity to determine thresholds"
325 " etc."
326 }, 1.e-8, 0, 1.e20}
327 );
328
329 slv_param_char(parameters,IPOPT_PARAM_MU_STRATEGY
330 ,(SlvParameterInitChar){{"mu_strategy"
331 ,"Update strategy for barrier parameter",6
332 ,"Determines which barrier parameter update strategy is to be used."
333 " 'monotone' is the monotone (Fiacco-McCormick) strategy;"
334 " 'adaptive' is the adaptive update strategy."
335 }, "monotone"}, (char *[]){
336 "monotone","adaptive",NULL
337 }
338 );
339
340 slv_param_bool(parameters,IPOPT_PARAM_SAFEEVAL
341 ,(SlvParameterInitBool){{"safeeval"
342 ,"Use safe evaluation?",1
343 ,"Use 'safe' function evaluation routines (TRUE) or allow ASCEND to "
344 "throw SIGFPE errors which will then halt integration (FALSE)."
345 }, FALSE}
346 );
347
348 slv_param_char(parameters,IPOPT_PARAM_DERIVATIVE_TEST
349 ,(SlvParameterInitChar){{"derivative_test"
350 ,"Use Derivative Checker?",1
351 ,"A finite-difference derivative checker is provided by IPOPT, which"
352 " will check Jacobian and gradient functions ('first-order') or"
353 " all first-order derivatives as well as the Hessian matrix"
354 " ('second-order'). The default is to perform no checks ('none')."
355 }, "none"}, (char *[]){
356 "none","first-order","second-order",NULL
357 }
358 );
359
360 slv_param_char(parameters,IPOPT_PARAM_HESS_APPROX
361 ,(SlvParameterInitChar){{"hessian_approximation"
362 ,"Hessian calculation method",1
363 ,"Use either an exact Hessian matrix based on symbolic derivatives"
364 " computed from the equations in the model ('exact'), or else use"
365 " a limited-memory quasi-Newton approximation ('limited-memory')."
366 " The default is 'limited-memory'."
367 }, "limited-memory"}, (char *[]){
368 "exact","limited-memory",NULL
369 }
370 );
371
372
373 asc_assert(parameters->num_parms==IPOPT_PARAMS);
374
375 return 1;
376 }
377
378 static void ipopt_get_parameters(slv_system_t server, SlvClientToken asys
379 , slv_parameters_t *parameters
380 ){
381 IpoptSystem *sys;
382 UNUSED_PARAMETER(server);
383
384 sys = SYS(asys);
385 //if(check_system(sys)) return;
386 mem_copy_cast(&(sys->p),parameters,sizeof(slv_parameters_t));
387 }
388
389
390 static void ipopt_set_parameters(slv_system_t server, SlvClientToken asys
391 ,slv_parameters_t *parameters
392 ){
393 IpoptSystem *sys;
394 UNUSED_PARAMETER(server);
395
396 sys = SYS(asys);
397 //if (check_system(sys)) return;
398 mem_copy_cast(parameters,&(sys->p),sizeof(slv_parameters_t));
399 }
400
401 /*------------------------------------------------------------------------------
402 EVALUATION AND RESULT HOOK FUNCTIONS
403 */
404
405 /**
406 update the model with new 'x' vector.
407 @return 0 on success.
408 */
409 int ipopt_update_model(IpoptSystem *sys, const double *x){
410 unsigned j;
411
412 CONSOLE_DEBUG("...");
413
414 asc_assert(sys);
415 asc_assert(sys->vlist);
416
417 /* FIXME do we need to update any other stuff? */
418 for(j = 0; j < sys->n; ++j){
419 var_set_value(sys->vlist[j], x[j]);
420 }
421
422 return 0;
423 }
424
425 /** Function to evaluate the objective function f(x).
426 @return 1 on success, 0 on failure
427
428 @param n (in), the number of variables in the problem (dimension of 'x').
429 @param x (in), the values for the primal variables, 'x' , at which 'f(x)' is to be evaluated.
430 @param new_x (in), false if any evaluation method was previously called with the same values in 'x', true otherwise.
431 @param obj_value (out) the value of the objective function ('f(x)').
432 */
433 Bool ipopt_eval_f(Index n, Number *x, Bool new_x, Number *obj_value, void *user_data){
434 IpoptSystem *sys;
435 sys = SYS(user_data);
436 int res;
437
438 CONSOLE_DEBUG("ipopt_eval_f");
439
440 asc_assert(n==sys->n);
441 asc_assert(sys->obj!=NULL);
442
443 if(new_x){
444 res = ipopt_update_model(sys,x);
445 if(res)return 0; /* fail model update */
446 }
447
448 sys->calc_ok = TRUE;
449
450 *obj_value = relman_eval(sys->obj,&(sys->calc_ok),SLV_PARAM_BOOL(&(sys->p),IPOPT_PARAM_SAFEEVAL));
451
452 return sys->calc_ok;
453 }
454
455 /**
456 @return 1 on success
457 */
458 Bool ipopt_eval_grad_f(Index n, Number* x, Bool new_x, Number* grad_f, void *user_data){
459 IpoptSystem *sys;
460 sys = SYS(user_data);
461 int j, res, len;
462 int count;
463 double *derivatives;
464 int *variables;
465 static var_filter_t vfilter = {
466 VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR
467 ,VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR | VAR_FIXED
468 };
469
470 CONSOLE_DEBUG("ipopt_eval_grad_f");
471
472 asc_assert(n==sys->n);
473 asc_assert(sys->obj);
474
475 if(new_x){
476 res = ipopt_update_model(sys,x);
477 if(res)return 0; /* fail model update */
478 }
479
480
481 /* evaluate grad_f(x) somehow */
482 for(j=0; j<n; ++j){
483 grad_f[j] = 0;
484 }
485
486 len = rel_n_incidences(sys->obj);
487 variables = ASC_NEW_ARRAY(int,len);
488 derivatives = ASC_NEW_ARRAY(double,len);
489
490 CONSOLE_DEBUG("allocated variables,derivatives");
491
492 relman_diff2(
493 sys->obj,&vfilter,derivatives,variables
494 , &count,SLV_PARAM_BOOL(&(sys->p),IPOPT_PARAM_SAFEEVAL)
495 );
496
497 for(j=0; j<len; ++j){
498 grad_f[variables[j]] = derivatives[j];
499 }
500
501 if(variables)ASC_FREE(variables);
502 if(derivatives)ASC_FREE(derivatives);
503
504 CONSOLE_DEBUG("done ipopt_eval_grad_f");
505 return 1; /* success, presumably */
506 }
507
508 Bool ipopt_eval_g(Index n, Number* x, Bool new_x, Index m, Number *g, void *user_data){
509 IpoptSystem *sys;
510 sys = SYS(user_data);
511 int i, res;
512 struct rel_relation *rel;
513 int calc_ok = 1;
514
515 CONSOLE_DEBUG("ipopt_eval_g (n=%d, m=%d)",sys->n, sys->m);
516
517 asc_assert(n==sys->n);
518 asc_assert(m==sys->m);
519
520 if(new_x){
521 res = ipopt_update_model(sys,x);
522 if(res)return 0; /* fail model update */
523 }
524
525 for(i=0;i<m;++i){
526 CONSOLE_DEBUG("rel %d: %s.",i,(sys->rlist[0] == sys->obj ? "OBJECTIVE" : "constraint"));
527 }
528
529 /** @todo constraint rels are all relations except the objective rel. do we need to sort the objective to the end? */
530 for(i=0; i<m; ++i){
531 rel = sys->rlist[i];
532 asc_assert(rel!=NULL);
533 g[i] = relman_eval(rel, &calc_ok,SLV_PARAM_BOOL(&(sys->p),IPOPT_PARAM_SAFEEVAL));
534 CONSOLE_DEBUG("g[%d] = %f",i,g[i]);
535 }
536
537 return calc_ok; /* fail: not yet implemented */
538 }
539
540 Bool ipopt_eval_jac_g(Index n, Number* x, Bool new_x, Index m
541 , Index nele_jac, Index* iRow, Index *jCol, Number* values
542 , void *user_data
543 ){
544 IpoptSystem *sys;
545 sys = SYS(user_data);
546 int i,res;
547
548 CONSOLE_DEBUG("ipopt_eval_jac_g");
549
550 if(!iRow || !jCol){
551 CONSOLE_DEBUG("sparsity structure requested");
552 }
553
554 asc_assert(sys!=NULL);
555 asc_assert(n==sys->n);
556 asc_assert(nele_jac==sys->nnzJ);
557 asc_assert(m==sys->m);
558
559 if(new_x){
560 res = ipopt_update_model(sys,x);
561 if(res)return 0; /* fail model update */
562 }
563
564 /* loop through the constraints */
565 for(i=0; i<m; ++i){
566 /* get derivatives for constraint i */
567 /* insert the derivatives into the matrix in row i, columns j */
568 }
569
570 CONSOLE_DEBUG("done ipopt_eval_jac_g");
571
572 return 0; /* fail: not yet implemented */
573 }
574
575 Bool ipopt_eval_h(Index n, Number* x, Bool new_x
576 , Number obj_factor, Index m, Number* lambda
577 , Bool new_lambda, Index nele_hess, Index* iRow
578 , Index* jCol, Number* values
579 , void *user_data
580 ){
581
582 CONSOLE_DEBUG("ipopt_eval_h");
583
584 if(iRow != NULL){
585 asc_assert(jCol !=NULL);
586 asc_assert(x==NULL); asc_assert(lambda==NULL); asc_assert(values==NULL);
587
588 /* identify the sparsity structure of the Hessian (note: only the lower-
589 left part is required by IPOPT , because the Hessian is symmetric) */
590 CONSOLE_DEBUG("Analysing of Hessian matrix sparsity structure not implemented");
591
592 /*
593 for(i=0; i<nvars; ++i){
594 for(j=i; j<nvars; ++j){
595 if(r2(sys->obj, sys->vlist[i], sys->vlist[j]) != 0){
596 iRow[nele_hess] = i;
597 jCol[nele_hess] = j;
598 nele_hess++
599 }
600 }
601 }
602 */
603
604 }else{
605 asc_assert(jCol==NULL);
606 asc_assert(x!=NULL); asc_assert(lambda!=NULL); asc_assert(values!=NULL);
607
608 /* evaluate the Hessian matrix */
609 CONSOLE_DEBUG("Evaluation of Hessian matrix not implemented");
610 }
611
612 return 0; /* fail: not yet implemented */
613 }
614
615 /*------------------------------------------------------------------------------
616 SOLVE ROUTINES
617 */
618
619 static int ipopt_presolve(slv_system_t server, SlvClientToken asys){
620 IpoptSystem *sys;
621 int max, i;
622 struct var_variable *var;
623
624 CONSOLE_DEBUG("PRESOLVE");
625
626 sys = SYS(asys);
627 ipopt_iteration_begins(sys);
628 //check_system(sys);
629
630 asc_assert(sys->vlist && sys->rlist);
631
632 /** @todo work out if matrix creation is not again needed */
633
634 /** @todo slv_sort_rels_and_vars(server,&(sys->m),&(sys->n)); */
635
636
637 /* count the number of optimisation variables */
638 sys->n = 0;
639 for(i = 0; i < sys->vtot; i++){
640 var = sys->vlist[i];
641 if(var_apply_filter(var,&(sys->vfilt))){
642 sys->n++;
643 }
644 }
645
646 /* set all relations as being 'unsatisfied' to start with... */
647 for(i=0; i < sys->rtot; ++i){
648 rel_set_satisfied(sys->rlist[i],FALSE);
649 }
650
651 sys->obj = slv_get_obj_relation(server); /*may have changed objective*/
652 if(!sys->obj){
653 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No objective function was specified");
654 return -3;
655 }
656 CONSOLE_DEBUG("got objective rel %p",sys->obj);
657 /* @todo check if old_obj == obj ? */
658
659 /* TODO are there cases where these should be different? */
660 sys->m = sys->vtot - 1; /* minus the objective relation */
661 CONSOLE_DEBUG("Numbers of constraints = %d",sys->m);
662
663 /** @todo we need to move the objective relation to the end of the list */
664
665 CONSOLE_DEBUG("got objective rel %p",sys->obj);
666
667 /* calculate nnz for hessian matrix @todo FIXME */
668
669 if(strcmp(SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_HESS_APPROX),"exact")==0){
670 /** @todo fix rtot to be 'm' instead */
671 sys->nnzH = relman_hessian_count(sys->rlist, sys->rtot, &(sys->vfilt), &(sys->rfilt), &max);
672 }else{
673 CONSOLE_DEBUG("Skipping relman_hessian_count as hessian method is not exact.");
674 sys->nnzH = sys->n * sys->m;
675 }
676
677 /* need to provide sparsity structure for hessian matrix? */
678
679 #if 0
680 /** @SEE http://www.coin-or.org/Ipopt/documentation/node37.html */
681 ipopt_eval_h(number_of_variables, NULL/*x at which to evaluate*/, TRUE /* new x */
682 , 1.0/*obj_factor*/, number_of_constraints, lambda/* values of the constraint multipliers */
683 , TRUE /* new lambda */, 0 /* number of nonzero elements in the Hessian */, Index* iRow
684 , Index* jCol, Number* values
685 , void *user_data
686 );
687 #endif
688
689
690 max = relman_obj_direction(sys->obj);
691 if(max==-1){
692 CONSOLE_DEBUG("this is a MINIMIZE problem");
693 }else{
694 CONSOLE_DEBUG("this is a MAXIMIZE problem");
695 }
696
697 CONSOLE_DEBUG("got %d constraints and %d vars in system", sys->m, sys->n);
698 /* calculate number of non-zeros in the Jacobian matrix for the constraint equations */
699
700 /* @todo make sure objective rel moved to end */
701 sys->nnzJ = relman_jacobian_count(sys->rlist, sys->m, &(sys->vfilt), &(sys->rfilt), &max);
702
703 CONSOLE_DEBUG("got %d non-zeros in constraint Jacobian", sys->nnzJ);
704
705 /* need to provide sparsity structure for jacobian? */
706
707
708
709 #if 0
710 if(sys->presolved > 0) { /* system has been presolved before */
711 if(!conopt_dof_changed(sys) /*no changes in fixed or included flags*/
712 && sys->p.partition == sys->J.old_partition
713 && sys->obj == sys->old_obj
714 ){
715 matrix_creation_needed = 0;
716 CONOPT_CONSOLE_DEBUG("YOU JUST AVOIDED MATRIX DESTRUCTION/CREATION");
717 }
718 }
719 #endif
720
721 #if 0
722 // check all this...
723
724 sys->presolved = 1; /* full presolve recognized here */
725 sys->resolve = 0; /* initialize resolve flag */
726
727 sys->J.old_partition = sys->p.partition;
728 sys->old_obj = sys->obj;
729
730 slv_sort_rels_and_vars(server,&(sys->con.m),&(sys->con.n));
731 CONOPT_CONSOLE_DEBUG("FOUND %d CONSTRAINTS AND %d VARS",sys->con.m,sys->con.n);
732 if (sys->obj != NULL) {
733 CONOPT_CONSOLE_DEBUG("ADDING OBJECT AS A ROW");
734 sys->con.m++; /* treat objective as a row */
735 }
736
737 cntvect = ASC_NEW_ARRAY(int,COIDEF_Size());
738 COIDEF_Ini(cntvect);
739 sys->con.cntvect = cntvect;
740 CONOPT_CONSOLE_DEBUG("NUMBER OF CONSTRAINTS = %d",sys->con.m);
741 COIDEF_NumVar(cntvect, &(sys->con.n));
742 COIDEF_NumCon(cntvect, &(sys->con.m));
743 sys->con.nz = num_jacobian_nonzeros(sys, &(sys->con.maxrow));
744 COIDEF_NumNZ(cntvect, &(sys->con.nz));
745 COIDEF_NumNlNz(cntvect, &(sys->con.nz));
746
747 sys->con.base = 0;
748 COIDEF_Base(cntvect,&(sys->con.base));
749 COIDEF_ErrLim(cntvect, &(DOMLIM));
750 COIDEF_ItLim(cntvect, &(ITER_LIMIT));
751
752 if(sys->obj!=NULL){
753 sys->con.optdir = relman_obj_direction(sys->obj);
754 sys->con.objcon = sys->con.m - 1; /* objective will be last row */
755 CONOPT_CONSOLE_DEBUG("SETTING OBJECTIVE CONSTRAINT TO BE %d",sys->con.objcon);
756 }else{
757 sys->con.optdir = 0;
758 sys->con.objcon = 0;
759 }
760 COIDEF_OptDir(cntvect, &(sys->con.optdir));
761 COIDEF_ObjCon(cntvect, &(sys->con.objcon));
762
763 temp = 0;
764 COIDEF_StdOut(cntvect, &temp);
765
766 int debugfv = 1;
767 COIDEF_DebugFV(cntvect, &debugfv);
768
769 destroy_vectors(sys);
770 destroy_matrices(sys);
771 create_matrices(server,sys);
772 create_vectors(sys);
773
774 sys->s.block.current_reordered_block = -2;
775 }
776
777 //...
778
779 if( matrix_creation_needed ) {
780 destroy_array(sys->s.cost);
781 sys->s.cost = create_zero_array(sys->s.costsize,struct slv_block_cost);
782 for( ind = 0; ind < sys->s.costsize; ++ind ) {
783 sys->s.cost[ind].reorder_method = -1;
784 }
785 } else {
786 reset_cost(sys->s.cost,sys->s.costsize);
787 }
788
789 #endif
790
791 /* Reset status */
792 sys->s.iteration = 0;
793 sys->s.cpu_elapsed = 0.0;
794 sys->s.converged = sys->s.diverged = sys->s.inconsistent = FALSE;
795 sys->s.block.previous_total_size = 0;
796 sys->s.costsize = 1+sys->s.block.number_of;
797
798
799 /* set to go to first unconverged block */
800 sys->s.block.current_block = -1;
801 sys->s.block.current_size = 0;
802 sys->s.calc_ok = TRUE;
803 sys->s.block.iteration = 0;
804 sys->obj_val = MAXDOUBLE/2000.0;
805
806 update_status(sys);
807
808 ipopt_iteration_ends(sys);
809
810 CONSOLE_DEBUG("Reset status");
811
812 /* sys->s.cost[sys->s.block.number_of].time=sys->s.cpu_elapsed; */
813
814 ERROR_REPORTER_HERE(ASC_USER_SUCCESS,"presolve completed");
815 return 0;
816 }
817
818
819 static int ipopt_solve(slv_system_t server, SlvClientToken asys){
820 IpoptSystem *sys;
821 UNUSED_PARAMETER(server);
822 enum ApplicationReturnStatus status;
823 int ret, i, j;
824 struct var_variable *var;
825
826 sys = SYS(asys);
827
828 double *x, *x_L, *x_U, *g_L, *g_U, *mult_x_L, *mult_x_U;
829
830 CONSOLE_DEBUG("SOLVING: sys->n = %d, sys->m = %d...",sys->n,sys->m);
831 asc_assert(sys->n!=-1);
832
833 /* set the number of variables and allocate space for the bounds */
834 x_L = ASC_NEW_ARRAY(Number,sys->n);
835 x_U = ASC_NEW_ARRAY(Number,sys->n);
836
837 CONSOLE_DEBUG("SETTING BOUNDS...");
838
839 /* @todo set the values for the variable bounds */
840 int jj = 0;
841 for(j = 0; j < sys->vtot; j++){
842 CONSOLE_DEBUG("j = %d, vtot = %d, vlist = %p",j,sys->vtot,sys->vlist);
843 var = sys->vlist[j];
844 if(var_apply_filter(var,&(sys->vfilt))){
845 CONSOLE_DEBUG("setting x_L[%d] = %e",jj,var_lower_bound(var));
846 assert(jj<sys->n);
847 x_L[jj] = var_lower_bound(var);
848 CONSOLE_DEBUG("setting x_U[%d] = %e",jj,var_upper_bound(var));
849 x_U[jj] = var_upper_bound(var);
850 jj++;
851 }
852 }
853
854 CONSOLE_DEBUG("jj = %d, sys->n = %d", jj, sys->n);
855 assert(jj==sys->n);
856
857 /** @todo set bounds on the constraints? */
858 /* is it possible to identify f(x)<a; f(x) >b and fold them into one? */
859 /* then find the constant parts and make then g_L or g_U accordingly */
860 /* what to do about other bounds? */
861 /* set the number of variables and allocate space for the bounds */
862 g_L = ASC_NEW_ARRAY(Number,sys->m);
863 g_U = ASC_NEW_ARRAY(Number,sys->m);
864 for(j = 0; j < sys->m; j++){
865 g_L[j] = 0;
866 g_U[j] = 0;
867 CONSOLE_DEBUG("set g_L[%d] = %e",j,g_L[j]);
868 CONSOLE_DEBUG("set g_U[%d] = %e",j,g_U[j]);
869 }
870
871 CONSOLE_DEBUG("CREATING PROBLEM...");
872
873 /* create the IpoptProblem */
874 CONSOLE_DEBUG("n = %d, m = %d, nnzJ = %d, nnzH = %d",sys->n, sys->m, sys->nnzJ, sys->nnzH);
875 sys->nlp = CreateIpoptProblem(sys->n, x_L, x_U, sys->m, g_L, g_U, sys->nnzJ, sys->nnzH, 0/*index style=C*/,
876 &ipopt_eval_f, &ipopt_eval_g, &ipopt_eval_grad_f,
877 &ipopt_eval_jac_g, &ipopt_eval_h
878 );
879
880 CONSOLE_DEBUG("FREEING INTERNAL STUFF");
881
882 /* We can free the memory now - the values for the bounds have been
883 copied internally in CreateIpoptProblem */
884 #if 0
885 /* freeing this stuff seems to cause a crash...?!?!? */
886 ASC_FREE(x_L);
887 ASC_FREE(x_U);
888 ASC_FREE(g_L);
889 ASC_FREE(g_U);
890 #endif
891
892 CONSOLE_DEBUG("SETTING OPTIONS...");
893 /* set some options */
894 AddIpoptNumOption(sys->nlp, "tol", SLV_PARAM_REAL(&(sys->p),IPOPT_PARAM_TOL));
895 AddIpoptStrOption(sys->nlp, "mu_strategy", SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_MU_STRATEGY));
896 AddIpoptStrOption(sys->nlp, "derivative_test", SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_DERIVATIVE_TEST));
897 AddIpoptStrOption(sys->nlp, "hessian_approximation", SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_HESS_APPROX));
898 AddIpoptStrOption(sys->nlp, "linear_solver", SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_LINEAR_SOLVER));
899
900 CONSOLE_DEBUG("Hessian method: %s",SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_HESS_APPROX));
901
902 /* initial values */
903 x = ASC_NEW_ARRAY(Number, sys->n);
904
905 /** @todo get values of 'x' from the model */
906
907 /* allocate space to store the bound multipliers at the solution */
908 mult_x_L = ASC_NEW_ARRAY(Number, sys->n);
909 mult_x_U = ASC_NEW_ARRAY(Number, sys->n);
910
911 CONSOLE_DEBUG("Calling IpoptSolve...");
912
913 /* solve the problem */
914 status = IpoptSolve(sys->nlp, x, NULL, &sys->obj_val, NULL, mult_x_L, mult_x_U, (void*)sys);
915
916 CONSOLE_DEBUG("Done IpoptSolve...");
917
918 /** @todo update the sys->s.xxxxx flags based on value of 'status' */
919
920 if (status == Solve_Succeeded) {
921 CONSOLE_DEBUG("Solution of the primal variables, x");
922 for (i=0; i<sys->n; i++) {
923 CONSOLE_DEBUG(" x[%d] = %e\n", i, x[i]);
924 }
925
926 CONSOLE_DEBUG("Solution of the bound multipliers, z_L and z_U");
927 for (i=0; i<sys->n; i++) {
928 CONSOLE_DEBUG(" z_L[%d] = %e", i, mult_x_L[i]);
929 }
930 for (i=0; i<sys->n; i++) {
931 CONSOLE_DEBUG(" z_U[%d] = %e", i, mult_x_U[i]);
932 }
933
934 CONSOLE_DEBUG("Objective value");
935 CONSOLE_DEBUG(" f(x*) = %e", sys->obj_val);
936
937 ret = 0; /* success */
938 }else{
939 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed solve, unknown status");
940 ret = 1; /* failure */
941 }
942
943 /* free allocated memory */
944 FreeIpoptProblem(sys->nlp);
945 ASC_FREE(x);
946 ASC_FREE(mult_x_L);
947 ASC_FREE(mult_x_U);
948
949 return ret;
950 }
951
952 /**
953 Prepare sys for entering an iteration, increasing the iteration counts
954 and starting the clock.
955 */
956 static void ipopt_iteration_begins(IpoptSystem *sys){
957 sys->clock = tm_cpu_time();
958 ++(sys->s.block.iteration);
959 ++(sys->s.iteration);
960 }
961
962
963 /*
964 Prepare sys for exiting an iteration, stopping the clock and recording
965 the cpu time.
966 */
967 static void ipopt_iteration_ends(IpoptSystem *sys){
968 double cpu_elapsed; /* elapsed this iteration */
969
970 cpu_elapsed = (double)(tm_cpu_time() - sys->clock);
971 sys->s.block.cpu_elapsed += cpu_elapsed;
972 sys->s.cpu_elapsed += cpu_elapsed;
973 }
974
975
976
977 static int ipopt_iterate(slv_system_t server, SlvClientToken asys){
978 CONSOLE_DEBUG("ipopt_iterate about to call ipopt_solve...");
979 return ipopt_solve(server,asys);
980 }
981
982 static int ipopt_resolve(slv_system_t server, SlvClientToken asys){
983 UNUSED_PARAMETER(server);
984
985 /* if implementing this, use the 'warm start' thing in IPOPT */
986
987 /* provide initial values of the 'multipliers' */
988
989 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
990 return 1;
991 }
992
993 static const SlvFunctionsT ipopt_internals = {
994 67
995 ,"IPOPT"
996 ,ipopt_create
997 ,ipopt_destroy
998 ,ipopt_eligible_solver
999 ,ipopt_get_default_parameters
1000 ,ipopt_get_parameters
1001 ,ipopt_set_parameters
1002 ,ipopt_get_status
1003 ,ipopt_solve
1004 ,ipopt_presolve
1005 ,ipopt_iterate
1006 ,ipopt_resolve
1007 ,NULL
1008 ,NULL
1009 ,NULL
1010 };
1011
1012 int ipopt_register(void){
1013 return solver_register(&ipopt_internals);
1014 }

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