/[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 1737 - (show annotations) (download) (as text)
Thu Feb 7 23:58:06 2008 UTC (16 years, 4 months ago) by jpye
File MIME type: text/x-csrc
File size: 27014 byte(s)
Further work on the IPOPT implementation.
1 /* ASCEND modelling environment
2 Copyright (C) 2007 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. I'VE ACTUALLY
23 ONLY JUST STARTED WRITING IT by starting with asc_tron.c and modifying.
24
25 The IPOPT solver is documented at http://projects.coin-or.org/Ipopt/
26 *//*
27 ASCEND wrapper for IPOPT originally by John Pye, Jun 2007.
28 */
29
30 #include <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 <solver/solver.h>
37
38 #include <system/calc.h>
39 #include <system/relman.h>
40 #include <system/slv_stdcalls.h>
41 #include <system/block.h>
42
43 #include <utilities/ascConfig.h>
44 #include <utilities/ascPanic.h>
45 #include <utilities/ascMalloc.h>
46 #include <utilities/ascDynaLoad.h>
47 #include <utilities/mem.h>
48 #include <utilities/ascEnvVar.h>
49 #include <general/tm_time.h>
50 #include <general/env.h>
51
52 #include <ipopt/IpStdCInterface.h>
53
54 ASC_DLLSPEC SolverRegisterFn ipopt_register;
55
56 /*------------------------------------------------------------------------------
57 DATA STRUCTURES AND FORWARD DEFS
58 */
59
60 /**
61 Documentation of solver options for IPOPT is at
62 http://www.coin-or.org/Ipopt/documentation/node1.html
63 */
64 enum{
65 IPOPT_PARAM_TOL
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 */
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 slv_param_real(parameters,IPOPT_PARAM_TOL
299 ,(SlvParameterInitReal){{"tol"
300 ,"Desired convergence tolerance (relative)",2
301 ,"Determines the convergence tolerance for the algorithm. The algorithm"
302 " terminates successfully, if the (scaled) NLP error becomes smaller"
303 " than this value, and if the (absolute) criteria according to "
304 " 'dual_inf_tol', 'primal_inf_tol', and 'cmpl_inf_tol' are met. (This"
305 " is epsilon_tol in Eqn. (6) in implementation paper). See also "
306 " 'acceptable_tol' as a second termination criterion. Note, some other"
307 " algorithmic features also use this quantity to determine thresholds"
308 " etc."
309 }, 1.e-8, 0, 1.e20}
310 );
311
312 slv_param_char(parameters,IPOPT_PARAM_MU_STRATEGY
313 ,(SlvParameterInitChar){{"mu_strategy"
314 ,"Update strategy for barrier parameter",6
315 ,"Determines which barrier parameter update strategy is to be used."
316 " 'monotone' is the monotone (Fiacco-McCormick) strategy;"
317 " 'adaptive' is the adaptive update strategy."
318 }, "monotone"}, (char *[]){
319 "monotone","adaptive",NULL
320 }
321 );
322
323 slv_param_bool(parameters,IPOPT_PARAM_SAFEEVAL
324 ,(SlvParameterInitBool){{"safeeval"
325 ,"Use safe evaluation?",1
326 ,"Use 'safe' function evaluation routines (TRUE) or allow ASCEND to "
327 "throw SIGFPE errors which will then halt integration (FALSE)."
328 }, FALSE}
329 );
330
331 slv_param_char(parameters,IPOPT_PARAM_DERIVATIVE_TEST
332 ,(SlvParameterInitChar){{"derivative_test"
333 ,"Use Derivative Checker?",1
334 ,"A finite-difference derivative checker is provided by IPOPT, which"
335 " will check Jacobian and gradient functions ('first-order') or"
336 " all first-order derivatives as well as the Hessian matrix"
337 " ('second-order'). The default is to perform no checks ('none')."
338 }, "none"}, (char *[]){
339 "none","first-order","second-order",NULL
340 }
341 );
342
343 slv_param_char(parameters,IPOPT_PARAM_HESS_APPROX
344 ,(SlvParameterInitChar){{"hessian_approximation"
345 ,"Hessian calculation method",1
346 ,"Use either an exact Hessian matrix based on symbolic derivatives"
347 " computed from the equations in the model ('exact'), or else use"
348 " a limited-memory quasi-Newton approximation ('limited-memory')."
349 " The default is 'exact'."
350 }, "exact"}, (char *[]){
351 "exact","limited-memory",NULL
352 }
353 );
354
355
356 asc_assert(parameters->num_parms==IPOPT_PARAMS);
357
358 return 1;
359 }
360
361 static void ipopt_get_parameters(slv_system_t server, SlvClientToken asys
362 , slv_parameters_t *parameters
363 ){
364 IpoptSystem *sys;
365 UNUSED_PARAMETER(server);
366
367 sys = SYS(asys);
368 //if(check_system(sys)) return;
369 mem_copy_cast(&(sys->p),parameters,sizeof(slv_parameters_t));
370 }
371
372
373 static void ipopt_set_parameters(slv_system_t server, SlvClientToken asys
374 ,slv_parameters_t *parameters
375 ){
376 IpoptSystem *sys;
377 UNUSED_PARAMETER(server);
378
379 sys = SYS(asys);
380 //if (check_system(sys)) return;
381 mem_copy_cast(parameters,&(sys->p),sizeof(slv_parameters_t));
382 }
383
384 /*------------------------------------------------------------------------------
385 EVALUATION AND RESULT HOOK FUNCTIONS
386 */
387
388 /**
389 update the model with new 'x' vector.
390 @return 0 on success.
391 */
392 int ipopt_update_model(IpoptSystem *sys, const double *x){
393 unsigned j;
394
395 asc_assert(sys);
396 asc_assert(sys->vlist);
397
398 /* FIXME do we need to update any other stuff? */
399 for(j = 0; j < sys->n; ++j){
400 var_set_value(sys->vlist[j], x[j]);
401 }
402
403 return 0;
404 }
405
406 /** Function to evaluate the objective function f(x).
407 @return 1 on success, 0 on failure
408
409 @param n (in), the number of variables in the problem (dimension of 'x').
410 @param x (in), the values for the primal variables, 'x' , at which 'f(x)' is to be evaluated.
411 @param new_x (in), false if any evaluation method was previously called with the same values in 'x', true otherwise.
412 @param obj_value (out) the value of the objective function ('f(x)').
413 */
414 Bool ipopt_eval_f(Index n, Number *x, Bool new_x, Number *obj_value, void *user_data){
415 IpoptSystem *sys;
416 sys = SYS(user_data);
417 int res;
418
419 CONSOLE_DEBUG("ipopt_eval_f");
420
421 asc_assert(n==sys->n);
422 asc_assert(sys->obj!=NULL);
423
424 if(new_x){
425 res = ipopt_update_model(sys,x);
426 if(res)return 0; /* fail model update */
427 }
428
429 sys->calc_ok = TRUE;
430
431 *obj_value = relman_eval(sys->obj,&(sys->calc_ok),SLV_PARAM_BOOL(&(sys->p),IPOPT_PARAM_SAFEEVAL));
432
433 return sys->calc_ok;
434 }
435
436 /**
437 @return 1 on success
438 */
439 Bool ipopt_eval_grad_f(Index n, Number* x, Bool new_x, Number* grad_f, void *user_data){
440 IpoptSystem *sys;
441 sys = SYS(user_data);
442 int j, res, len;
443 int count;
444 double *derivatives;
445 int *variables;
446 static var_filter_t vfilter = {
447 VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR
448 ,VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR | VAR_FIXED
449 };
450
451 CONSOLE_DEBUG("ipopt_eval_grad_f");
452
453 asc_assert(n==sys->n);
454 asc_assert(sys->obj);
455
456 if(new_x){
457 res = ipopt_update_model(sys,x);
458 if(res)return 0; /* fail model update */
459 }
460
461
462 /* evaluate grad_f(x) somehow */
463 for(j=0; j<n; ++j){
464 grad_f[j] = 0;
465 }
466
467 len = rel_n_incidences(sys->obj);
468 variables = ASC_NEW_ARRAY(int,len);
469 derivatives = ASC_NEW_ARRAY(double,len);
470
471 relman_diff2(
472 sys->obj,&vfilter,derivatives,variables
473 , &count,SLV_PARAM_BOOL(&(sys->p),IPOPT_PARAM_SAFEEVAL)
474 );
475
476 for(j=0; j<len; ++j){
477 grad_f[variables[j]] = derivatives[j];
478 }
479
480 ASC_FREE(variables);
481 ASC_FREE(derivatives);
482
483 return 1; /* success, presumably */
484 }
485
486 Bool ipopt_eval_g(Index n, Number* x, Bool new_x, Index m, Number *g, void *user_data){
487 IpoptSystem *sys;
488 sys = SYS(user_data);
489 int i, res;
490
491 CONSOLE_DEBUG("ipopt_eval_g");
492
493 asc_assert(n==sys->n);
494 asc_assert(m==sys->m);
495
496 if(new_x){
497 res = ipopt_update_model(sys,x);
498 if(res)return 0; /* fail model update */
499 }
500
501 for(i=0; i<m; ++i){
502 g[i] = 0;
503 }
504
505 return 0; /* fail: not yet implemented */
506 }
507
508 Bool ipopt_eval_jac_g(Index n, Number* x, Bool new_x, Index m
509 , Index nele_jac, Index* iRow, Index *jCol, Number* values
510 , void *user_data
511 ){
512 IpoptSystem *sys;
513 sys = SYS(user_data);
514 int i,res;
515
516 asc_assert(sys!=NULL);
517 asc_assert(n==sys->n);
518 asc_assert(nele_jac==sys->nnzJ);
519 asc_assert(m==sys->m);
520
521 if(new_x){
522 res = ipopt_update_model(sys,x);
523 if(res)return 0; /* fail model update */
524 }
525
526 /* loop through the constraints */
527 for(i=0; i<m; ++i){
528 /* get derivatives for constraint i */
529 /* insert the derivatives into the matrix in row i, columns j */
530 }
531
532 return 0; /* fail: not yet implemented */
533 }
534
535 Bool ipopt_eval_h(Index n, Number* x, Bool new_x
536 , Number obj_factor, Index m, Number* lambda
537 , Bool new_lambda, Index nele_hess, Index* iRow
538 , Index* jCol, Number* values
539 , void *user_data
540 ){
541
542 CONSOLE_DEBUG("ipopt_eval_h");
543
544 if(iRow != NULL){
545 asc_assert(jCol !=NULL);
546 asc_assert(x==NULL); asc_assert(lambda==NULL); asc_assert(values==NULL);
547
548 /* identify the sparsity structure of the Hessian (note: only the lower-
549 left part is required by IPOPT , because the Hessian is symmetric) */
550 CONSOLE_DEBUG("Analysing of Hessian matrix sparsity structure not implemented");
551
552 /*
553 for(i=0; i<nvars; ++i){
554 for(j=i; j<nvars; ++j){
555 if(r2(sys->obj, sys->vlist[i], sys->vlist[j]) != 0){
556 iRow[nele_hess] = i;
557 jCol[nele_hess] = j;
558 nele_hess++
559 }
560 }
561 }
562 */
563
564 }else{
565 asc_assert(jCol==NULL);
566 asc_assert(x!=NULL); asc_assert(lambda!=NULL); asc_assert(values!=NULL);
567
568 /* evaluate the Hessian matrix */
569 CONSOLE_DEBUG("Evaluation of Hessian matrix not implemented");
570 }
571
572 return 0; /* fail: not yet implemented */
573 }
574
575 /*------------------------------------------------------------------------------
576 SOLVE ROUTINES
577 */
578
579 static int ipopt_presolve(slv_system_t server, SlvClientToken asys){
580 IpoptSystem *sys;
581 int max, i;
582
583 CONSOLE_DEBUG("PRESOLVE");
584
585 sys = SYS(asys);
586 ipopt_iteration_begins(sys);
587 //check_system(sys);
588
589 asc_assert(sys->vlist && sys->rlist);
590
591 /** @TODO work out if matrix creation is not again needed */
592
593 /** @TODO slv_sort_rels_and_vars(server,&(sys->m),&(sys->n)); */
594
595 /* TODO are there cases where these should be different? */
596 sys->n = sys->rtot;
597 sys->m = sys->vtot;
598
599 /* set all relations as being 'unsatisfied' to start with... */
600 for(i=0; i < sys->rtot; ++i){
601 rel_set_satisfied(sys->rlist[i],FALSE);
602 }
603
604 sys->obj = slv_get_obj_relation(server); /*may have changed objective*/
605 if(!sys->obj){
606 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No objective function was specified");
607 return -3;
608 }
609
610 CONSOLE_DEBUG("got objective rel %p",sys->obj);
611
612 /* calculate nnz for hessian matrix @TODO FIXME */
613
614 if(strcmp(SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_HESS_APPROX),"exact")==0){
615 sys->nnzH = relman_hessian_count(sys->rlist, sys->rtot, &(sys->vfilt), &(sys->rfilt), &max);
616 }else{
617 CONSOLE_DEBUG("Skipping relman_hessian_count as hessian method is not exact.");
618 sys->nnzH = sys->n * sys->m;
619 }
620
621 /* need to provide sparsity structure for hessian matrix? */
622
623 #if 0
624 /** @SEE http://www.coin-or.org/Ipopt/documentation/node37.html */
625 ipopt_eval_h(number_of_variables, NULL/*x at which to evaluate*/, TRUE /* new x */
626 , 1.0/*obj_factor*/, number_of_constraints, lambda/* values of the constraint multipliers */
627 , TRUE /* new lambda */, 0 /* number of nonzero elements in the Hessian */, Index* iRow
628 , Index* jCol, Number* values
629 , void *user_data
630 );
631 #endif
632
633
634 max = relman_obj_direction(sys->obj);
635 if(max==-1){
636 CONSOLE_DEBUG("this is a MINIMIZE problem");
637 }else{
638 CONSOLE_DEBUG("this is a MAXIMIZE problem");
639 }
640
641 CONSOLE_DEBUG("got %d relations and %d vars in system", sys->rtot, sys->vtot);
642 /* calculate number of non-zeros in the Jacobian matrix for the constraint equations */
643
644 sys->nnzJ = relman_jacobian_count(sys->rlist, sys->rtot, &(sys->vfilt), &(sys->rfilt), &max);
645
646 CONSOLE_DEBUG("got %d non-zeros in constraint Jacobian", sys->nnzJ);
647
648 /* need to provide sparsity structure for jacobian? */
649
650
651
652 #if 0
653 if(sys->presolved > 0) { /* system has been presolved before */
654 if(!conopt_dof_changed(sys) /*no changes in fixed or included flags*/
655 && sys->p.partition == sys->J.old_partition
656 && sys->obj == sys->old_obj
657 ){
658 matrix_creation_needed = 0;
659 CONOPT_CONSOLE_DEBUG("YOU JUST AVOIDED MATRIX DESTRUCTION/CREATION");
660 }
661 }
662 #endif
663
664 #if 0
665 // check all this...
666
667 sys->presolved = 1; /* full presolve recognized here */
668 sys->resolve = 0; /* initialize resolve flag */
669
670 sys->J.old_partition = sys->p.partition;
671 sys->old_obj = sys->obj;
672
673 slv_sort_rels_and_vars(server,&(sys->con.m),&(sys->con.n));
674 CONOPT_CONSOLE_DEBUG("FOUND %d CONSTRAINTS AND %d VARS",sys->con.m,sys->con.n);
675 if (sys->obj != NULL) {
676 CONOPT_CONSOLE_DEBUG("ADDING OBJECT AS A ROW");
677 sys->con.m++; /* treat objective as a row */
678 }
679
680 cntvect = ASC_NEW_ARRAY(int,COIDEF_Size());
681 COIDEF_Ini(cntvect);
682 sys->con.cntvect = cntvect;
683 CONOPT_CONSOLE_DEBUG("NUMBER OF CONSTRAINTS = %d",sys->con.m);
684 COIDEF_NumVar(cntvect, &(sys->con.n));
685 COIDEF_NumCon(cntvect, &(sys->con.m));
686 sys->con.nz = num_jacobian_nonzeros(sys, &(sys->con.maxrow));
687 COIDEF_NumNZ(cntvect, &(sys->con.nz));
688 COIDEF_NumNlNz(cntvect, &(sys->con.nz));
689
690 sys->con.base = 0;
691 COIDEF_Base(cntvect,&(sys->con.base));
692 COIDEF_ErrLim(cntvect, &(DOMLIM));
693 COIDEF_ItLim(cntvect, &(ITER_LIMIT));
694
695 if(sys->obj!=NULL){
696 sys->con.optdir = relman_obj_direction(sys->obj);
697 sys->con.objcon = sys->con.m - 1; /* objective will be last row */
698 CONOPT_CONSOLE_DEBUG("SETTING OBJECTIVE CONSTRAINT TO BE %d",sys->con.objcon);
699 }else{
700 sys->con.optdir = 0;
701 sys->con.objcon = 0;
702 }
703 COIDEF_OptDir(cntvect, &(sys->con.optdir));
704 COIDEF_ObjCon(cntvect, &(sys->con.objcon));
705
706 temp = 0;
707 COIDEF_StdOut(cntvect, &temp);
708
709 int debugfv = 1;
710 COIDEF_DebugFV(cntvect, &debugfv);
711
712 destroy_vectors(sys);
713 destroy_matrices(sys);
714 create_matrices(server,sys);
715 create_vectors(sys);
716
717 sys->s.block.current_reordered_block = -2;
718 }
719
720 //...
721
722 if( matrix_creation_needed ) {
723 destroy_array(sys->s.cost);
724 sys->s.cost = create_zero_array(sys->s.costsize,struct slv_block_cost);
725 for( ind = 0; ind < sys->s.costsize; ++ind ) {
726 sys->s.cost[ind].reorder_method = -1;
727 }
728 } else {
729 reset_cost(sys->s.cost,sys->s.costsize);
730 }
731
732 #endif
733
734 /* Reset status */
735 sys->s.iteration = 0;
736 sys->s.cpu_elapsed = 0.0;
737 sys->s.converged = sys->s.diverged = sys->s.inconsistent = FALSE;
738 sys->s.block.previous_total_size = 0;
739 sys->s.costsize = 1+sys->s.block.number_of;
740
741
742 /* set to go to first unconverged block */
743 sys->s.block.current_block = -1;
744 sys->s.block.current_size = 0;
745 sys->s.calc_ok = TRUE;
746 sys->s.block.iteration = 0;
747 sys->obj_val = MAXDOUBLE/2000.0;
748
749 update_status(sys);
750
751 ipopt_iteration_ends(sys);
752
753 CONSOLE_DEBUG("Reset status");
754
755 /* sys->s.cost[sys->s.block.number_of].time=sys->s.cpu_elapsed; */
756
757 ERROR_REPORTER_HERE(ASC_USER_SUCCESS,"presolve completed");
758 return 0;
759 }
760
761
762 static int ipopt_solve(slv_system_t server, SlvClientToken asys){
763 IpoptSystem *sys;
764 UNUSED_PARAMETER(server);
765 enum ApplicationReturnStatus status;
766 int ret, i, j;
767 struct var_variable *var;
768
769 sys = SYS(asys);
770
771 double *x, *x_L, *x_U, *g_L, *g_U, *mult_x_L, *mult_x_U;
772
773 CONSOLE_DEBUG("SOLVING, sys->n = %d...",sys->n);
774 asc_assert(sys->n!=-1);
775
776 /* set the number of variables and allocate space for the bounds */
777 x_L = ASC_NEW_ARRAY(Number,sys->n);
778 x_U = ASC_NEW_ARRAY(Number,sys->n);
779
780 CONSOLE_DEBUG("SETTING BOUNDS...");
781
782 /* @TODO set the values for the variable bounds */
783 int jj = 0;
784 for(j = 0; j < sys->vtot; j++){
785 CONSOLE_DEBUG("j = %d, vtot = %d, vlist = %p",j,sys->vtot,sys->vlist);
786 var = sys->vlist[j];
787 if(var_apply_filter(var,&(sys->vfilt))){
788 CONSOLE_DEBUG("setting x_L[%d]",jj);
789 x_L[jj] = var_lower_bound(var);
790 x_U[jj] = var_upper_bound(var);
791 jj++;
792 }
793 }
794
795 /** @TODO set bounds on the constraints? */
796 /* need to identify equations that share the same non-constant parts? */
797 /* then find the constant parts and make then g_L or g_U accordingly */
798 /* what to do about other bounds? */
799 /* set the number of variables and allocate space for the bounds */
800 g_L = ASC_NEW_ARRAY(Number,sys->m);
801 g_U = ASC_NEW_ARRAY(Number,sys->m);
802 for(j = 0; j < sys->m; j++){
803 g_L[j] = 0;
804 g_U[j] = 0;
805 }
806
807 CONSOLE_DEBUG("CREATING PROBLEM...");
808
809 /* create the IpoptProblem */
810 CONSOLE_DEBUG("n = %d, m = %d, nnzJ = %d, nnzH = %d",sys->n, sys->m, sys->nnzJ, sys->nnzH);
811 sys->nlp = CreateIpoptProblem(sys->n, x_L, x_U, sys->m, g_L, g_U, sys->nnzJ, sys->nnzH, 0/*index style=C*/,
812 &ipopt_eval_f, &ipopt_eval_g, &ipopt_eval_grad_f,
813 &ipopt_eval_jac_g, &ipopt_eval_h
814 );
815
816 CONSOLE_DEBUG("FREEING INTERNAL STUFF");
817
818 /* We can free the memory now - the values for the bounds have been
819 copied internally in CreateIpoptProblem */
820 #if 0
821 /* freeing this stuff seems to cause a crash...?!?!? */
822 ASC_FREE(x_L);
823 ASC_FREE(x_U);
824 ASC_FREE(g_L);
825 ASC_FREE(g_U);
826 #endif
827
828 CONSOLE_DEBUG("SETTING OPTIONS...");
829 /* set some options */
830 AddIpoptNumOption(sys->nlp, "tol", SLV_PARAM_REAL(&(sys->p),IPOPT_PARAM_TOL));
831 AddIpoptStrOption(sys->nlp, "mu_strategy", SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_MU_STRATEGY));
832 AddIpoptStrOption(sys->nlp, "derivative_test", SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_DERIVATIVE_TEST));
833 AddIpoptStrOption(sys->nlp, "hessian_approximation", SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_HESS_APPROX));
834
835 CONSOLE_DEBUG("Hessian method: %s",SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_HESS_APPROX));
836
837 /* initial values */
838 x = ASC_NEW_ARRAY(Number, sys->n);
839
840 /** @TODO get values of 'x' from the model */
841
842 /* allocate space to store the bound multipliers at the solution */
843 mult_x_L = ASC_NEW_ARRAY(Number, sys->n);
844 mult_x_U = ASC_NEW_ARRAY(Number, sys->n);
845
846 CONSOLE_DEBUG("Calling IpoptSolve...");
847
848 /* solve the problem */
849 status = IpoptSolve(sys->nlp, x, NULL, &sys->obj_val, NULL, mult_x_L, mult_x_U, (void*)sys);
850
851 CONSOLE_DEBUG("Done IpoptSolve...");
852
853 /** @TODO update the sys->s.xxxxx flags based on value of 'status' */
854
855 if (status == Solve_Succeeded) {
856 CONSOLE_DEBUG("Solution of the primal variables, x");
857 for (i=0; i<sys->n; i++) {
858 CONSOLE_DEBUG(" x[%d] = %e\n", i, x[i]);
859 }
860
861 CONSOLE_DEBUG("Solution of the bound multipliers, z_L and z_U");
862 for (i=0; i<sys->n; i++) {
863 CONSOLE_DEBUG(" z_L[%d] = %e", i, mult_x_L[i]);
864 }
865 for (i=0; i<sys->n; i++) {
866 CONSOLE_DEBUG(" z_U[%d] = %e", i, mult_x_U[i]);
867 }
868
869 CONSOLE_DEBUG("Objective value");
870 CONSOLE_DEBUG(" f(x*) = %e", sys->obj_val);
871
872 ret = 0; /* success */
873 }else{
874 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed solve, unknown status");
875 ret = 1; /* failure */
876 }
877
878 /* free allocated memory */
879 FreeIpoptProblem(sys->nlp);
880 ASC_FREE(x);
881 ASC_FREE(mult_x_L);
882 ASC_FREE(mult_x_U);
883
884 return ret;
885 }
886
887 /**
888 Prepare sys for entering an iteration, increasing the iteration counts
889 and starting the clock.
890 */
891 static void ipopt_iteration_begins(IpoptSystem *sys){
892 sys->clock = tm_cpu_time();
893 ++(sys->s.block.iteration);
894 ++(sys->s.iteration);
895 }
896
897
898 /*
899 Prepare sys for exiting an iteration, stopping the clock and recording
900 the cpu time.
901 */
902 static void ipopt_iteration_ends(IpoptSystem *sys){
903 double cpu_elapsed; /* elapsed this iteration */
904
905 cpu_elapsed = (double)(tm_cpu_time() - sys->clock);
906 sys->s.block.cpu_elapsed += cpu_elapsed;
907 sys->s.cpu_elapsed += cpu_elapsed;
908 }
909
910
911
912 static int ipopt_iterate(slv_system_t server, SlvClientToken asys){
913 CONSOLE_DEBUG("ipopt_iterate about to call ipopt_solve...");
914 return ipopt_solve(server,asys);
915 }
916
917 static int ipopt_resolve(slv_system_t server, SlvClientToken asys){
918 UNUSED_PARAMETER(server);
919
920 /* if implementing this, use the 'warm start' thing in IPOPT */
921
922 /* provide initial values of the 'multipliers' */
923
924 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
925 return 1;
926 }
927
928 static const SlvFunctionsT ipopt_internals = {
929 67
930 ,"IPOPT"
931 ,ipopt_create
932 ,ipopt_destroy
933 ,ipopt_eligible_solver
934 ,ipopt_get_default_parameters
935 ,ipopt_get_parameters
936 ,ipopt_set_parameters
937 ,ipopt_get_status
938 ,ipopt_solve
939 ,ipopt_presolve
940 ,ipopt_iterate
941 ,ipopt_resolve
942 ,NULL
943 ,NULL
944 ,NULL
945 };
946
947 int ipopt_register(void){
948 return solver_register(&ipopt_internals);
949 }

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