/[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 1771 - (show annotations) (download) (as text)
Mon May 12 04:28:04 2008 UTC (12 years, 4 months ago) by jpye
File MIME type: text/x-csrc
File size: 27289 byte(s)
Debug IPOPT. Seems that next thing to implement eval_jac_g sparsity structure.
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 <coin/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 CONSOLE_DEBUG("allocated variables,derivatives");
472
473 relman_diff2(
474 sys->obj,&vfilter,derivatives,variables
475 , &count,SLV_PARAM_BOOL(&(sys->p),IPOPT_PARAM_SAFEEVAL)
476 );
477
478 for(j=0; j<len; ++j){
479 grad_f[variables[j]] = derivatives[j];
480 }
481
482 if(variables)ASC_FREE(variables);
483 if(derivatives)ASC_FREE(derivatives);
484
485 CONSOLE_DEBUG("done ipopt_eval_grad_f");
486 return 1; /* success, presumably */
487 }
488
489 Bool ipopt_eval_g(Index n, Number* x, Bool new_x, Index m, Number *g, void *user_data){
490 IpoptSystem *sys;
491 sys = SYS(user_data);
492 int i, res;
493
494 CONSOLE_DEBUG("ipopt_eval_g");
495
496 asc_assert(n==sys->n);
497 asc_assert(m==sys->m);
498
499 if(new_x){
500 res = ipopt_update_model(sys,x);
501 if(res)return 0; /* fail model update */
502 }
503
504 for(i=0; i<m; ++i){
505 g[i] = 0;
506 }
507
508 return 0; /* fail: not yet implemented */
509 }
510
511 Bool ipopt_eval_jac_g(Index n, Number* x, Bool new_x, Index m
512 , Index nele_jac, Index* iRow, Index *jCol, Number* values
513 , void *user_data
514 ){
515 IpoptSystem *sys;
516 sys = SYS(user_data);
517 int i,res;
518
519 CONSOLE_DEBUG("ipopt_eval_jac_g");
520
521 if(!iRow || !jCol){
522 CONSOLE_DEBUG("sparsity structure requested");
523 }
524
525 asc_assert(sys!=NULL);
526 asc_assert(n==sys->n);
527 asc_assert(nele_jac==sys->nnzJ);
528 asc_assert(m==sys->m);
529
530 if(new_x){
531 res = ipopt_update_model(sys,x);
532 if(res)return 0; /* fail model update */
533 }
534
535 /* loop through the constraints */
536 for(i=0; i<m; ++i){
537 /* get derivatives for constraint i */
538 /* insert the derivatives into the matrix in row i, columns j */
539 }
540
541 CONSOLE_DEBUG("done ipopt_eval_jac_g");
542
543 return 0; /* fail: not yet implemented */
544 }
545
546 Bool ipopt_eval_h(Index n, Number* x, Bool new_x
547 , Number obj_factor, Index m, Number* lambda
548 , Bool new_lambda, Index nele_hess, Index* iRow
549 , Index* jCol, Number* values
550 , void *user_data
551 ){
552
553 CONSOLE_DEBUG("ipopt_eval_h");
554
555 if(iRow != NULL){
556 asc_assert(jCol !=NULL);
557 asc_assert(x==NULL); asc_assert(lambda==NULL); asc_assert(values==NULL);
558
559 /* identify the sparsity structure of the Hessian (note: only the lower-
560 left part is required by IPOPT , because the Hessian is symmetric) */
561 CONSOLE_DEBUG("Analysing of Hessian matrix sparsity structure not implemented");
562
563 /*
564 for(i=0; i<nvars; ++i){
565 for(j=i; j<nvars; ++j){
566 if(r2(sys->obj, sys->vlist[i], sys->vlist[j]) != 0){
567 iRow[nele_hess] = i;
568 jCol[nele_hess] = j;
569 nele_hess++
570 }
571 }
572 }
573 */
574
575 }else{
576 asc_assert(jCol==NULL);
577 asc_assert(x!=NULL); asc_assert(lambda!=NULL); asc_assert(values!=NULL);
578
579 /* evaluate the Hessian matrix */
580 CONSOLE_DEBUG("Evaluation of Hessian matrix not implemented");
581 }
582
583 return 0; /* fail: not yet implemented */
584 }
585
586 /*------------------------------------------------------------------------------
587 SOLVE ROUTINES
588 */
589
590 static int ipopt_presolve(slv_system_t server, SlvClientToken asys){
591 IpoptSystem *sys;
592 int max, i;
593
594 CONSOLE_DEBUG("PRESOLVE");
595
596 sys = SYS(asys);
597 ipopt_iteration_begins(sys);
598 //check_system(sys);
599
600 asc_assert(sys->vlist && sys->rlist);
601
602 /** @TODO work out if matrix creation is not again needed */
603
604 /** @TODO slv_sort_rels_and_vars(server,&(sys->m),&(sys->n)); */
605
606 /* TODO are there cases where these should be different? */
607 sys->n = sys->rtot;
608 sys->m = sys->vtot;
609
610 /* set all relations as being 'unsatisfied' to start with... */
611 for(i=0; i < sys->rtot; ++i){
612 rel_set_satisfied(sys->rlist[i],FALSE);
613 }
614
615 sys->obj = slv_get_obj_relation(server); /*may have changed objective*/
616 if(!sys->obj){
617 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No objective function was specified");
618 return -3;
619 }
620
621 CONSOLE_DEBUG("got objective rel %p",sys->obj);
622
623 /* calculate nnz for hessian matrix @TODO FIXME */
624
625 if(strcmp(SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_HESS_APPROX),"exact")==0){
626 sys->nnzH = relman_hessian_count(sys->rlist, sys->rtot, &(sys->vfilt), &(sys->rfilt), &max);
627 }else{
628 CONSOLE_DEBUG("Skipping relman_hessian_count as hessian method is not exact.");
629 sys->nnzH = sys->n * sys->m;
630 }
631
632 /* need to provide sparsity structure for hessian matrix? */
633
634 #if 0
635 /** @SEE http://www.coin-or.org/Ipopt/documentation/node37.html */
636 ipopt_eval_h(number_of_variables, NULL/*x at which to evaluate*/, TRUE /* new x */
637 , 1.0/*obj_factor*/, number_of_constraints, lambda/* values of the constraint multipliers */
638 , TRUE /* new lambda */, 0 /* number of nonzero elements in the Hessian */, Index* iRow
639 , Index* jCol, Number* values
640 , void *user_data
641 );
642 #endif
643
644
645 max = relman_obj_direction(sys->obj);
646 if(max==-1){
647 CONSOLE_DEBUG("this is a MINIMIZE problem");
648 }else{
649 CONSOLE_DEBUG("this is a MAXIMIZE problem");
650 }
651
652 CONSOLE_DEBUG("got %d relations and %d vars in system", sys->rtot, sys->vtot);
653 /* calculate number of non-zeros in the Jacobian matrix for the constraint equations */
654
655 sys->nnzJ = relman_jacobian_count(sys->rlist, sys->rtot, &(sys->vfilt), &(sys->rfilt), &max);
656
657 CONSOLE_DEBUG("got %d non-zeros in constraint Jacobian", sys->nnzJ);
658
659 /* need to provide sparsity structure for jacobian? */
660
661
662
663 #if 0
664 if(sys->presolved > 0) { /* system has been presolved before */
665 if(!conopt_dof_changed(sys) /*no changes in fixed or included flags*/
666 && sys->p.partition == sys->J.old_partition
667 && sys->obj == sys->old_obj
668 ){
669 matrix_creation_needed = 0;
670 CONOPT_CONSOLE_DEBUG("YOU JUST AVOIDED MATRIX DESTRUCTION/CREATION");
671 }
672 }
673 #endif
674
675 #if 0
676 // check all this...
677
678 sys->presolved = 1; /* full presolve recognized here */
679 sys->resolve = 0; /* initialize resolve flag */
680
681 sys->J.old_partition = sys->p.partition;
682 sys->old_obj = sys->obj;
683
684 slv_sort_rels_and_vars(server,&(sys->con.m),&(sys->con.n));
685 CONOPT_CONSOLE_DEBUG("FOUND %d CONSTRAINTS AND %d VARS",sys->con.m,sys->con.n);
686 if (sys->obj != NULL) {
687 CONOPT_CONSOLE_DEBUG("ADDING OBJECT AS A ROW");
688 sys->con.m++; /* treat objective as a row */
689 }
690
691 cntvect = ASC_NEW_ARRAY(int,COIDEF_Size());
692 COIDEF_Ini(cntvect);
693 sys->con.cntvect = cntvect;
694 CONOPT_CONSOLE_DEBUG("NUMBER OF CONSTRAINTS = %d",sys->con.m);
695 COIDEF_NumVar(cntvect, &(sys->con.n));
696 COIDEF_NumCon(cntvect, &(sys->con.m));
697 sys->con.nz = num_jacobian_nonzeros(sys, &(sys->con.maxrow));
698 COIDEF_NumNZ(cntvect, &(sys->con.nz));
699 COIDEF_NumNlNz(cntvect, &(sys->con.nz));
700
701 sys->con.base = 0;
702 COIDEF_Base(cntvect,&(sys->con.base));
703 COIDEF_ErrLim(cntvect, &(DOMLIM));
704 COIDEF_ItLim(cntvect, &(ITER_LIMIT));
705
706 if(sys->obj!=NULL){
707 sys->con.optdir = relman_obj_direction(sys->obj);
708 sys->con.objcon = sys->con.m - 1; /* objective will be last row */
709 CONOPT_CONSOLE_DEBUG("SETTING OBJECTIVE CONSTRAINT TO BE %d",sys->con.objcon);
710 }else{
711 sys->con.optdir = 0;
712 sys->con.objcon = 0;
713 }
714 COIDEF_OptDir(cntvect, &(sys->con.optdir));
715 COIDEF_ObjCon(cntvect, &(sys->con.objcon));
716
717 temp = 0;
718 COIDEF_StdOut(cntvect, &temp);
719
720 int debugfv = 1;
721 COIDEF_DebugFV(cntvect, &debugfv);
722
723 destroy_vectors(sys);
724 destroy_matrices(sys);
725 create_matrices(server,sys);
726 create_vectors(sys);
727
728 sys->s.block.current_reordered_block = -2;
729 }
730
731 //...
732
733 if( matrix_creation_needed ) {
734 destroy_array(sys->s.cost);
735 sys->s.cost = create_zero_array(sys->s.costsize,struct slv_block_cost);
736 for( ind = 0; ind < sys->s.costsize; ++ind ) {
737 sys->s.cost[ind].reorder_method = -1;
738 }
739 } else {
740 reset_cost(sys->s.cost,sys->s.costsize);
741 }
742
743 #endif
744
745 /* Reset status */
746 sys->s.iteration = 0;
747 sys->s.cpu_elapsed = 0.0;
748 sys->s.converged = sys->s.diverged = sys->s.inconsistent = FALSE;
749 sys->s.block.previous_total_size = 0;
750 sys->s.costsize = 1+sys->s.block.number_of;
751
752
753 /* set to go to first unconverged block */
754 sys->s.block.current_block = -1;
755 sys->s.block.current_size = 0;
756 sys->s.calc_ok = TRUE;
757 sys->s.block.iteration = 0;
758 sys->obj_val = MAXDOUBLE/2000.0;
759
760 update_status(sys);
761
762 ipopt_iteration_ends(sys);
763
764 CONSOLE_DEBUG("Reset status");
765
766 /* sys->s.cost[sys->s.block.number_of].time=sys->s.cpu_elapsed; */
767
768 ERROR_REPORTER_HERE(ASC_USER_SUCCESS,"presolve completed");
769 return 0;
770 }
771
772
773 static int ipopt_solve(slv_system_t server, SlvClientToken asys){
774 IpoptSystem *sys;
775 UNUSED_PARAMETER(server);
776 enum ApplicationReturnStatus status;
777 int ret, i, j;
778 struct var_variable *var;
779
780 sys = SYS(asys);
781
782 double *x, *x_L, *x_U, *g_L, *g_U, *mult_x_L, *mult_x_U;
783
784 CONSOLE_DEBUG("SOLVING, sys->n = %d...",sys->n);
785 asc_assert(sys->n!=-1);
786
787 /* set the number of variables and allocate space for the bounds */
788 x_L = ASC_NEW_ARRAY(Number,sys->n);
789 x_U = ASC_NEW_ARRAY(Number,sys->n);
790
791 CONSOLE_DEBUG("SETTING BOUNDS...");
792
793 /* @TODO set the values for the variable bounds */
794 int jj = 0;
795 for(j = 0; j < sys->vtot; j++){
796 CONSOLE_DEBUG("j = %d, vtot = %d, vlist = %p",j,sys->vtot,sys->vlist);
797 var = sys->vlist[j];
798 if(var_apply_filter(var,&(sys->vfilt))){
799 CONSOLE_DEBUG("setting x_L[%d]",jj);
800 x_L[jj] = var_lower_bound(var);
801 x_U[jj] = var_upper_bound(var);
802 jj++;
803 }
804 }
805
806 /** @TODO set bounds on the constraints? */
807 /* need to identify equations that share the same non-constant parts? */
808 /* then find the constant parts and make then g_L or g_U accordingly */
809 /* what to do about other bounds? */
810 /* set the number of variables and allocate space for the bounds */
811 g_L = ASC_NEW_ARRAY(Number,sys->m);
812 g_U = ASC_NEW_ARRAY(Number,sys->m);
813 for(j = 0; j < sys->m; j++){
814 g_L[j] = 0;
815 g_U[j] = 0;
816 }
817
818 CONSOLE_DEBUG("CREATING PROBLEM...");
819
820 /* create the IpoptProblem */
821 CONSOLE_DEBUG("n = %d, m = %d, nnzJ = %d, nnzH = %d",sys->n, sys->m, sys->nnzJ, sys->nnzH);
822 sys->nlp = CreateIpoptProblem(sys->n, x_L, x_U, sys->m, g_L, g_U, sys->nnzJ, sys->nnzH, 0/*index style=C*/,
823 &ipopt_eval_f, &ipopt_eval_g, &ipopt_eval_grad_f,
824 &ipopt_eval_jac_g, &ipopt_eval_h
825 );
826
827 CONSOLE_DEBUG("FREEING INTERNAL STUFF");
828
829 /* We can free the memory now - the values for the bounds have been
830 copied internally in CreateIpoptProblem */
831 #if 0
832 /* freeing this stuff seems to cause a crash...?!?!? */
833 ASC_FREE(x_L);
834 ASC_FREE(x_U);
835 ASC_FREE(g_L);
836 ASC_FREE(g_U);
837 #endif
838
839 CONSOLE_DEBUG("SETTING OPTIONS...");
840 /* set some options */
841 AddIpoptNumOption(sys->nlp, "tol", SLV_PARAM_REAL(&(sys->p),IPOPT_PARAM_TOL));
842 AddIpoptStrOption(sys->nlp, "mu_strategy", SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_MU_STRATEGY));
843 AddIpoptStrOption(sys->nlp, "derivative_test", SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_DERIVATIVE_TEST));
844 AddIpoptStrOption(sys->nlp, "hessian_approximation", SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_HESS_APPROX));
845
846 CONSOLE_DEBUG("Hessian method: %s",SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_HESS_APPROX));
847
848 /* initial values */
849 x = ASC_NEW_ARRAY(Number, sys->n);
850
851 /** @TODO get values of 'x' from the model */
852
853 /* allocate space to store the bound multipliers at the solution */
854 mult_x_L = ASC_NEW_ARRAY(Number, sys->n);
855 mult_x_U = ASC_NEW_ARRAY(Number, sys->n);
856
857 CONSOLE_DEBUG("Calling IpoptSolve...");
858
859 /* solve the problem */
860 status = IpoptSolve(sys->nlp, x, NULL, &sys->obj_val, NULL, mult_x_L, mult_x_U, (void*)sys);
861
862 CONSOLE_DEBUG("Done IpoptSolve...");
863
864 /** @TODO update the sys->s.xxxxx flags based on value of 'status' */
865
866 if (status == Solve_Succeeded) {
867 CONSOLE_DEBUG("Solution of the primal variables, x");
868 for (i=0; i<sys->n; i++) {
869 CONSOLE_DEBUG(" x[%d] = %e\n", i, x[i]);
870 }
871
872 CONSOLE_DEBUG("Solution of the bound multipliers, z_L and z_U");
873 for (i=0; i<sys->n; i++) {
874 CONSOLE_DEBUG(" z_L[%d] = %e", i, mult_x_L[i]);
875 }
876 for (i=0; i<sys->n; i++) {
877 CONSOLE_DEBUG(" z_U[%d] = %e", i, mult_x_U[i]);
878 }
879
880 CONSOLE_DEBUG("Objective value");
881 CONSOLE_DEBUG(" f(x*) = %e", sys->obj_val);
882
883 ret = 0; /* success */
884 }else{
885 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed solve, unknown status");
886 ret = 1; /* failure */
887 }
888
889 /* free allocated memory */
890 FreeIpoptProblem(sys->nlp);
891 ASC_FREE(x);
892 ASC_FREE(mult_x_L);
893 ASC_FREE(mult_x_U);
894
895 return ret;
896 }
897
898 /**
899 Prepare sys for entering an iteration, increasing the iteration counts
900 and starting the clock.
901 */
902 static void ipopt_iteration_begins(IpoptSystem *sys){
903 sys->clock = tm_cpu_time();
904 ++(sys->s.block.iteration);
905 ++(sys->s.iteration);
906 }
907
908
909 /*
910 Prepare sys for exiting an iteration, stopping the clock and recording
911 the cpu time.
912 */
913 static void ipopt_iteration_ends(IpoptSystem *sys){
914 double cpu_elapsed; /* elapsed this iteration */
915
916 cpu_elapsed = (double)(tm_cpu_time() - sys->clock);
917 sys->s.block.cpu_elapsed += cpu_elapsed;
918 sys->s.cpu_elapsed += cpu_elapsed;
919 }
920
921
922
923 static int ipopt_iterate(slv_system_t server, SlvClientToken asys){
924 CONSOLE_DEBUG("ipopt_iterate about to call ipopt_solve...");
925 return ipopt_solve(server,asys);
926 }
927
928 static int ipopt_resolve(slv_system_t server, SlvClientToken asys){
929 UNUSED_PARAMETER(server);
930
931 /* if implementing this, use the 'warm start' thing in IPOPT */
932
933 /* provide initial values of the 'multipliers' */
934
935 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
936 return 1;
937 }
938
939 static const SlvFunctionsT ipopt_internals = {
940 67
941 ,"IPOPT"
942 ,ipopt_create
943 ,ipopt_destroy
944 ,ipopt_eligible_solver
945 ,ipopt_get_default_parameters
946 ,ipopt_get_parameters
947 ,ipopt_set_parameters
948 ,ipopt_get_status
949 ,ipopt_solve
950 ,ipopt_presolve
951 ,ipopt_iterate
952 ,ipopt_resolve
953 ,NULL
954 ,NULL
955 ,NULL
956 };
957
958 int ipopt_register(void){
959 return solver_register(&ipopt_internals);
960 }

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