/[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 1555 - (show annotations) (download) (as text)
Tue Jul 24 14:16:45 2007 UTC (13 years, 3 months ago) by jpye
File MIME type: text/x-csrc
File size: 22942 byte(s)
Bit more work on IPOPT presolve and hessian calculation.
Added export of relman_jacobian_count.
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_PARAMS
70 };
71
72 #define SYS(s) ((IpoptSystem *)(s))
73
74 struct IpoptSystemStruct{
75
76 /*
77 Problem definition
78 */
79 slv_system_t slv; /* slv_system_t back-link */
80 struct rel_relation *obj; /* Objective function: NULL = none */
81 struct rel_relation *old_obj;/* Objective function: NULL = none */
82 struct var_variable **vlist; /* Variable list (NULL terminated) */
83 struct rel_relation **rlist; /* Relation list (NULL terminated) */
84
85 var_filter_t vfilt;
86 rel_filter_t rfilt;
87
88 /*
89 Solver information
90 */
91 int32 presolved; /* ? Has the system been presolved */
92 int32 resolve; /* ? Has the system been resolved */
93 slv_parameters_t p; /* Parameters */
94 slv_status_t s; /* Status (as of iteration end) */
95
96 int32 cap; /* Order of matrix/vectors */
97 int32 rank; /* Symbolic rank of problem */
98 int32 vused; /* Free and incident variables */
99 int32 vtot; /* length of varlist */
100 int32 rused; /* Included relations */
101 int32 rtot; /* length of rellist */
102 double clock; /* CPU time */
103
104 int32 calc_ok;
105 double obj_val;
106
107 void *parm_array[IPOPT_PARAMS];
108 struct slv_parameter pa[IPOPT_PARAMS];
109
110 /*
111 IPOPT DATA
112 */
113 Index n; /* number of variables */
114 Index m; /* number of constraints */
115
116 int nnzJ; /* number of non zeros in the jacobian of the constraints */
117 int nnzH; /* number of non-zeros in the hessian of the objective */
118
119 Number* x_L; /* lower bounds on x */
120 Number* x_U; /* upper bounds on x */
121 Number* g_L; /* lower bounds on g */
122 Number* g_U; /* upper bounds on g */
123
124 IpoptProblem nlp; /* IpoptProblem */
125
126 enum ApplicationReturnStatus status; /* Solve return code */
127 Number* x; /* starting point and solution vector */
128 Number* mult_x_L; /* lower bound multipliers at the solution */
129 Number* mult_x_U; /* upper bound multipliers at the solution */
130 Index i; /* generic counter */
131 };
132
133 typedef struct IpoptSystemStruct IpoptSystem;
134
135 static int ipopt_get_default_parameters(slv_system_t server, SlvClientToken asys, slv_parameters_t *parameters);
136
137 static void ipopt_iteration_begins(IpoptSystem *sys);
138 static void ipopt_iteration_ends(IpoptSystem *sys);
139
140 /*------------------------------------------------------------------------------
141 SYSTEM SETUP/DESTROY, STATUS AND SOLVER ELIGIBILITY
142 */
143
144 static SlvClientToken ipopt_create(slv_system_t server, int32*statusindex){
145 IpoptSystem *sys;
146
147 sys = ASC_NEW_CLEAR(IpoptSystem);
148 if(sys==NULL){
149 *statusindex = 1;
150 return sys;
151 }
152
153 sys->p.parms = sys->pa;
154 sys->p.dynamic_parms = 0;
155 ipopt_get_default_parameters(server,(SlvClientToken)sys,&(sys->p));
156 sys->p.whose = (*statusindex);
157
158 sys->presolved = 0;
159 sys->resolve = 0;
160
161 sys->n = -1;
162 sys->m = -1;
163
164 sys->s.ok = TRUE;
165 sys->s.calc_ok = TRUE;
166 sys->s.costsize = 0;
167 sys->s.cost = NULL; /*redundant, but sanity preserving */
168 sys->s.block.number_of = 1;
169 sys->s.block.current_block = 0;
170 sys->s.block.current_reordered_block = 0;
171 sys->s.block.current_size = 0;
172 sys->s.block.previous_total_size = 0;
173 sys->s.block.iteration = 0;
174 sys->s.block.funcs = 0;
175 sys->s.block.jacs = 0;
176 sys->s.block.cpu_elapsed = 0;
177 sys->s.block.functime = 0;
178 sys->s.block.jactime = 0;
179 sys->s.block.residual = 0;
180
181 sys->rfilt.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
182 sys->rfilt.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
183 sys->vfilt.matchbits = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR | VAR_FIXED);
184 sys->vfilt.matchvalue = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR);
185
186 sys->vlist = slv_get_solvers_var_list(server);
187 sys->rlist = slv_get_solvers_rel_list(server);
188
189 sys->rtot = slv_get_num_solvers_rels(server);
190 sys->vtot = slv_get_num_solvers_vars(server);
191
192 sys->obj = slv_get_obj_relation(server);
193 if(sys->vlist == NULL) {
194 ASC_FREE(sys);
195 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IPOPT called with no variables.");
196 *statusindex = -2;
197 return NULL;
198 }
199 if(sys->rlist == NULL && sys->obj == NULL) {
200 ASC_FREE(sys);
201 ERROR_REPORTER_HERE(ASC_PROG_ERR,"IPOPT called with no relations or objective.");
202 *statusindex = -1;
203 return NULL;
204 }
205
206 /* do nothing with the objective list, pars, bounds, extrels, etc etc */
207
208 slv_check_var_initialization(server);
209 *statusindex = 0;
210 return((SlvClientToken)sys);
211 }
212
213 static int32 ipopt_destroy(slv_system_t server, SlvClientToken asys){
214 UNUSED_PARAMETER(server);
215 ERROR_REPORTER_HERE(ASC_PROG_ERR,"ipopt_destroy not implemented");
216 return 1;
217 }
218
219
220 static int ipopt_get_status(slv_system_t server, SlvClientToken asys
221 ,slv_status_t *status
222 ){
223 IpoptSystem *sys;
224 (void)server; /* stop gcc whine about unused parameter */
225
226 sys = SYS(asys);
227 //if (check_system(sys)) return 1;
228 mem_copy_cast(&(sys->s),status,sizeof(slv_status_t));
229 return 0;
230 }
231
232 static int32 ipopt_eligible_solver(slv_system_t server){
233 UNUSED_PARAMETER(server);
234
235 /// TODO check that there is a MAXIMIZE or MINIMIZE statement
236 /// TODO check that there are no discrete-valued variables
237 /// TODO check that there are no WHENs or CONDITIONALs
238 /// TODO check anything else?
239
240 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"ipopt_eligible_solver not implemented");
241 return 1;
242 }
243
244 /*------------------------------------------------------------------------------
245 SOLVER PARAMETERS
246 */
247
248 static
249 int32 ipopt_get_default_parameters(slv_system_t server, SlvClientToken asys
250 ,slv_parameters_t *parameters
251 ){
252 IpoptSystem *sys = NULL;
253 struct slv_parameter *new_parms = NULL;
254 int32 make_macros = 0;
255
256 if(server != NULL && asys != NULL) {
257 sys = SYS(asys);
258 make_macros = 1;
259 }
260
261 if(parameters->parms == NULL) {
262 new_parms = ASC_NEW_ARRAY_OR_NULL(struct slv_parameter,IPOPT_PARAMS);
263 if(new_parms == NULL) {
264 return -1;
265 }
266 parameters->parms = new_parms;
267 parameters->dynamic_parms = 1;
268 }
269
270 parameters->num_parms = 0;
271
272 slv_param_int(parameters,IPOPT_PARAM_MAX_ITER
273 ,(SlvParameterInitInt){{"max_iter"
274 ,"Maximum number of iterations",2
275 ,"The algorithm terminates with an error message if the number of iterations exceeded this number."
276 }, 3000, 0, 100000000}
277 );
278
279 slv_param_real(parameters,IPOPT_PARAM_TOL
280 ,(SlvParameterInitReal){{"tol"
281 ,"Desired convergence tolerance (relative)",2
282 ,"Determines the convergence tolerance for the algorithm. The algorithm"
283 " terminates successfully, if the (scaled) NLP error becomes smaller"
284 " than this value, and if the (absolute) criteria according to "
285 " 'dual_inf_tol', 'primal_inf_tol', and 'cmpl_inf_tol' are met. (This"
286 " is epsilon_tol in Eqn. (6) in implementation paper). See also "
287 " 'acceptable_tol' as a second termination criterion. Note, some other"
288 " algorithmic features also use this quantity to determine thresholds"
289 " etc."
290 }, 1e-8, 0, 1e20}
291 );
292
293 slv_param_char(parameters,IPOPT_PARAM_MU_STRATEGY
294 ,(SlvParameterInitChar){{"mu_strategy"
295 ,"Update strategy for barrier parameter",6
296 ,"Determines which barrier parameter update strategy is to be used."
297 " 'monotone' is the monotone (Fiacco-McCormick) strategy;"
298 " 'adaptive' is the adaptive update strategy."
299 }, "monotone"}, (char *[]){
300 "monotone","adaptive",NULL
301 }
302 );
303
304 slv_param_bool(parameters,IPOPT_PARAM_SAFEEVAL
305 ,(SlvParameterInitBool){{"safeeval"
306 ,"Use safe evaluation?",1
307 ,"Use 'safe' function evaluation routines (TRUE) or allow ASCEND to "
308 "throw SIGFPE errors which will then halt integration (FALSE)."
309 }, FALSE}
310 );
311
312 asc_assert(parameters->num_parms==IPOPT_PARAMS);
313
314 return 1;
315 }
316
317 static void ipopt_get_parameters(slv_system_t server, SlvClientToken asys
318 , slv_parameters_t *parameters
319 ){
320 IpoptSystem *sys;
321 UNUSED_PARAMETER(server);
322
323 sys = SYS(asys);
324 //if(check_system(sys)) return;
325 mem_copy_cast(&(sys->p),parameters,sizeof(slv_parameters_t));
326 }
327
328
329 static void ipopt_set_parameters(slv_system_t server, SlvClientToken asys
330 ,slv_parameters_t *parameters
331 ){
332 IpoptSystem *sys;
333 UNUSED_PARAMETER(server);
334
335 sys = SYS(asys);
336 //if (check_system(sys)) return;
337 mem_copy_cast(parameters,&(sys->p),sizeof(slv_parameters_t));
338 }
339
340 /*------------------------------------------------------------------------------
341 EVALUATION AND RESULT HOOK FUNCTIONS
342 */
343
344 /**
345 update the model with new 'x' vector.
346 @return 0 on success.
347 */
348 int ipopt_update_model(IpoptSystem *sys, const double *x){
349 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
350 return 1; /* error */
351 }
352
353 /** Function to evaluate the objective function f(x).
354 @return 1 on success, 0 on failure
355
356 @param n (in), the number of variables in the problem (dimension of 'x').
357 @param x (in), the values for the primal variables, 'x' , at which 'f(x)' is to be evaluated.
358 @param new_x (in), false if any evaluation method was previously called with the same values in 'x', true otherwise.
359 @param obj_value (out) the value of the objective function ('f(x)').
360 */
361 Bool ipopt_eval_f(Index n, Number *x, Bool new_x, Number *obj_value, void *user_data){
362 IpoptSystem *sys;
363 sys = SYS(user_data);
364 int res;
365
366 asc_assert(n==sys->n);
367 asc_assert(sys->obj!=NULL);
368
369 if(new_x){
370 res = ipopt_update_model(sys,x);
371 if(res)return 0; /* fail model update */
372 }
373
374 sys->calc_ok = TRUE;
375
376 *obj_value = relman_eval(sys->obj,&(sys->calc_ok),SLV_PARAM_BOOL(&(sys->p),IPOPT_PARAM_SAFEEVAL));
377
378 return sys->calc_ok;
379 }
380
381 /**
382 @return 1 on success
383 */
384 Bool ipopt_eval_grad_f(Index n, Number* x, Bool new_x, Number* grad_f, void *user_data){
385 IpoptSystem *sys;
386 sys = SYS(user_data);
387 int j, res, len;
388 int count;
389 double *derivatives;
390 int *variables;
391 static var_filter_t vfilter = {
392 VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR
393 ,VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR | VAR_FIXED
394 };
395
396 asc_assert(n==sys->n);
397 asc_assert(sys->obj);
398
399 if(new_x){
400 res = ipopt_update_model(sys,x);
401 if(res)return 0; /* fail model update */
402 }
403
404
405 /* evaluate grad_f(x) somehow */
406 for(j=0; j<n; ++j){
407 grad_f[j] = 0;
408 }
409
410 len = rel_n_incidences(sys->obj);
411 variables = ASC_NEW_ARRAY(int,len);
412 derivatives = ASC_NEW_ARRAY(double,len);
413
414 relman_diff2(
415 sys->obj,&vfilter,derivatives,variables
416 , &count,SLV_PARAM_BOOL(&(sys->p),IPOPT_PARAM_SAFEEVAL)
417 );
418
419 for(j=0; j<len; ++j){
420 grad_f[variables[j]] = derivatives[j];
421 }
422
423 ASC_FREE(variables);
424 ASC_FREE(derivatives);
425
426 return 1; /* success, presumably */
427 }
428
429 Bool ipopt_eval_g(Index n, Number* x, Bool new_x, Index m, Number *g, void *user_data){
430 IpoptSystem *sys;
431 sys = SYS(user_data);
432 int i, res;
433
434 asc_assert(n==sys->n);
435 asc_assert(m==sys->m);
436
437 if(new_x){
438 res = ipopt_update_model(sys,x);
439 if(res)return 0; /* fail model update */
440 }
441
442 for(i=0; i<m; ++i){
443 g[i] = 0;
444 }
445
446 return 0; /* fail: not yet implemented */
447 }
448
449 Bool ipopt_eval_jac_g(Index n, Number* x, Bool new_x, Index m
450 , Index nele_jac, Index* iRow, Index *jCol, Number* values
451 , void *user_data
452 ){
453 IpoptSystem *sys;
454 sys = SYS(user_data);
455 int i,res;
456
457 asc_assert(n==sys->n);
458 asc_assert(nele_jac==sys->nnzJ);
459 asc_assert(m==sys->m);
460
461 if(new_x){
462 res = ipopt_update_model(sys,x);
463 if(res)return 0; /* fail model update */
464 }
465
466 /* loop through the constraints */
467 for(i=0; i<m; ++i){
468 /* get derivatives for constraint i */
469 /* insert the derivatives into the matrix in row i, columns j */
470 }
471
472 return 0; /* fail: not yet implemented */
473 }
474
475 Bool ipopt_eval_h(Index n, Number* x, Bool new_x
476 , Number obj_factor, Index m, Number* lambda
477 , Bool new_lambda, Index nele_hess, Index* iRow
478 , Index* jCol, Number* values
479 , void *user_data
480 ){
481 if(iRow != NULL){
482 asc_assert(jCol !=NULL);
483 asc_assert(x==NULL); asc_assert(lambda==NULL); asc_assert(values==NULL);
484
485 /* identify the sparsity structure of the Hessian (note: only the lower-
486 left part is required by IPOPT , because the Hessian is symmetric) */
487
488 /*
489 for(i=0; i<nvars; ++i){
490 for(j=i; j<nvars; ++j){
491 if(relman_hess_expr(sys->obj, sys->vlist[i], sys->vlist[j]) != 0){
492 iRow[nele_hess] = i;
493 jCol[nele_hess] = j;
494 nele_hess++
495 }
496 }
497 }
498 */
499
500 }else{
501 asc_assert(jCol==NULL);
502 asc_assert(x!=NULL); asc_assert(lambda!=NULL); asc_assert(values!=NULL);
503
504 /* evaluate the Hessian matrix */
505 }
506
507 return 0; /* fail: not yet implemented */
508 }
509
510 /*------------------------------------------------------------------------------
511 SOLVE ROUTINES
512 */
513
514 static int ipopt_presolve(slv_system_t server, SlvClientToken asys){
515 IpoptSystem *sys;
516 int max, i;
517
518 CONSOLE_DEBUG("PRESOLVE");
519
520 sys = SYS(asys);
521 ipopt_iteration_begins(sys);
522 //check_system(sys);
523
524 asc_assert(sys->vlist && sys->rlist);
525
526 /** @TODO work out if matrix creation is not again needed */
527
528 /** @TODO slv_sort_rels_and_vars(server,&(sys->m),&(sys->n)); */
529
530 /* set all relations as being 'unsatisfied' to start with... */
531 for(i=0; i < sys->rtot; ++i){
532 rel_set_satisfied(sys->rlist[i],FALSE);
533 }
534
535 sys->obj = slv_get_obj_relation(server); /*may have changed objective*/
536 if(!sys->obj){
537 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No objective function was specified");
538 return -3;
539 }
540
541 CONSOLE_DEBUG("got objective rel %p",sys->obj);
542
543 /** @TODO calculate nnz for hessian matrix */
544
545 /* need to provide sparsity structure for hessian matrix? */
546
547 max = relman_obj_direction(sys->obj);
548 if(max==-1){
549 CONSOLE_DEBUG("this is a MINIMIZE problem");
550 }else{
551 CONSOLE_DEBUG("this is a MAXIMIZE problem");
552 }
553
554 CONSOLE_DEBUG("got %d relations and %d vars in system", sys->rtot, sys->vtot);
555 /* calculate number of non-zeros in the Jacobian matrix for the constraint equations */
556
557 sys->nnzJ = relman_jacobian_count(sys->rlist, sys->rtot, &(sys->vfilt), &(sys->rfilt), &max);
558
559 CONSOLE_DEBUG("got %d non-zeros in constraint Jacobian", sys->nnzJ);
560
561 /* need to provide sparsity structure for jacobian? */
562
563
564
565 #if 0
566 if(sys->presolved > 0) { /* system has been presolved before */
567 if(!conopt_dof_changed(sys) /*no changes in fixed or included flags*/
568 && sys->p.partition == sys->J.old_partition
569 && sys->obj == sys->old_obj
570 ){
571 matrix_creation_needed = 0;
572 CONOPT_CONSOLE_DEBUG("YOU JUST AVOIDED MATRIX DESTRUCTION/CREATION");
573 }
574 }
575 #endif
576
577 #if 0
578 // check all this...
579
580 sys->presolved = 1; /* full presolve recognized here */
581 sys->resolve = 0; /* initialize resolve flag */
582
583 sys->J.old_partition = sys->p.partition;
584 sys->old_obj = sys->obj;
585
586 slv_sort_rels_and_vars(server,&(sys->con.m),&(sys->con.n));
587 CONOPT_CONSOLE_DEBUG("FOUND %d CONSTRAINTS AND %d VARS",sys->con.m,sys->con.n);
588 if (sys->obj != NULL) {
589 CONOPT_CONSOLE_DEBUG("ADDING OBJECT AS A ROW");
590 sys->con.m++; /* treat objective as a row */
591 }
592
593 cntvect = ASC_NEW_ARRAY(int,COIDEF_Size());
594 COIDEF_Ini(cntvect);
595 sys->con.cntvect = cntvect;
596 CONOPT_CONSOLE_DEBUG("NUMBER OF CONSTRAINTS = %d",sys->con.m);
597 COIDEF_NumVar(cntvect, &(sys->con.n));
598 COIDEF_NumCon(cntvect, &(sys->con.m));
599 sys->con.nz = num_jacobian_nonzeros(sys, &(sys->con.maxrow));
600 COIDEF_NumNZ(cntvect, &(sys->con.nz));
601 COIDEF_NumNlNz(cntvect, &(sys->con.nz));
602
603 sys->con.base = 0;
604 COIDEF_Base(cntvect,&(sys->con.base));
605 COIDEF_ErrLim(cntvect, &(DOMLIM));
606 COIDEF_ItLim(cntvect, &(ITER_LIMIT));
607
608 if(sys->obj!=NULL){
609 sys->con.optdir = relman_obj_direction(sys->obj);
610 sys->con.objcon = sys->con.m - 1; /* objective will be last row */
611 CONOPT_CONSOLE_DEBUG("SETTING OBJECTIVE CONSTRAINT TO BE %d",sys->con.objcon);
612 }else{
613 sys->con.optdir = 0;
614 sys->con.objcon = 0;
615 }
616 COIDEF_OptDir(cntvect, &(sys->con.optdir));
617 COIDEF_ObjCon(cntvect, &(sys->con.objcon));
618
619 temp = 0;
620 COIDEF_StdOut(cntvect, &temp);
621
622 int debugfv = 1;
623 COIDEF_DebugFV(cntvect, &debugfv);
624
625 destroy_vectors(sys);
626 destroy_matrices(sys);
627 create_matrices(server,sys);
628 create_vectors(sys);
629
630 sys->s.block.current_reordered_block = -2;
631 }
632
633 //...
634
635 if( matrix_creation_needed ) {
636 destroy_array(sys->s.cost);
637 sys->s.cost = create_zero_array(sys->s.costsize,struct slv_block_cost);
638 for( ind = 0; ind < sys->s.costsize; ++ind ) {
639 sys->s.cost[ind].reorder_method = -1;
640 }
641 } else {
642 reset_cost(sys->s.cost,sys->s.costsize);
643 }
644
645 #endif
646
647 /* Reset status */
648 sys->s.iteration = 0;
649 sys->s.cpu_elapsed = 0.0;
650 sys->s.converged = sys->s.diverged = sys->s.inconsistent = FALSE;
651 sys->s.block.previous_total_size = 0;
652 sys->s.costsize = 1+sys->s.block.number_of;
653
654
655 /* set to go to first unconverged block */
656 sys->s.block.current_block = -1;
657 sys->s.block.current_size = 0;
658 sys->s.calc_ok = TRUE;
659 sys->s.block.iteration = 0;
660 sys->obj_val = MAXDOUBLE/2000.0;
661
662 ipopt_iteration_ends(sys);
663
664 CONSOLE_DEBUG("Reset status");
665
666 /* sys->s.cost[sys->s.block.number_of].time=sys->s.cpu_elapsed; */
667
668 ERROR_REPORTER_HERE(ASC_PROG_ERR,"presolve completed");
669 return 0;
670 }
671
672
673 static int ipopt_solve(slv_system_t server, SlvClientToken asys){
674 IpoptSystem *sys;
675 UNUSED_PARAMETER(server);
676 enum ApplicationReturnStatus status;
677 int ret, i, j;
678 struct var_variable *var;
679
680 sys = SYS(asys);
681
682 double *x, *x_L, *x_U, *g_L, *g_U, *mult_x_L, *mult_x_U;
683
684 CONSOLE_DEBUG("SOLVING...");
685
686 /* set the number of variables and allocate space for the bounds */
687 x_L = ASC_NEW_ARRAY(Number,sys->n);
688 x_U = ASC_NEW_ARRAY(Number,sys->n);
689
690 /* @TODO set the values for the variable bounds */
691 int jj = 0;
692 for(j = 0; j < sys->vtot; j++){
693 var = sys->vlist[j];
694 if(var_apply_filter(var,&(sys->vfilt))){
695 x_L[jj] = var_lower_bound(var);
696 x_U[jj] = var_upper_bound(var);
697 jj++;
698 }
699 }
700
701 /** @TODO set bounds on the constraints? */
702
703 /* create the IpoptProblem */
704 sys->nlp = CreateIpoptProblem(sys->n, x_L, x_U, sys->m, g_L, g_U, sys->nnzJ, sys->nnzH, 0/*index style=C*/,
705 &ipopt_eval_f, &ipopt_eval_g, &ipopt_eval_grad_f,
706 &ipopt_eval_jac_g, &ipopt_eval_h
707 );
708
709 /* We can free the memory now - the values for the bounds have been
710 copied internally in CreateIpoptProblem */
711 ASC_FREE(x_L);
712 ASC_FREE(x_U);
713 ASC_FREE(g_L);
714 ASC_FREE(g_U);
715
716 /* set some options */
717 AddIpoptNumOption(sys->nlp, "tol", SLV_PARAM_BOOL(&(sys->p),IPOPT_PARAM_TOL));
718 AddIpoptStrOption(sys->nlp, "mu_strategy", SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_MU_STRATEGY));
719
720 /* initial values */
721 x = ASC_NEW_ARRAY(Number, sys->n);
722
723 /** @TODO get values of 'x' from the model */
724
725 /* allocate space to store the bound multipliers at the solution */
726 mult_x_L = ASC_NEW_ARRAY(Number, sys->n);
727 mult_x_U = ASC_NEW_ARRAY(Number, sys->n);
728
729 /* solve the problem */
730 status = IpoptSolve(sys->nlp, x, NULL, &sys->obj_val, NULL, mult_x_L, mult_x_U, NULL);
731
732 /** @TODO update the sys->s.xxxxx flags based on value of 'status' */
733
734 if (status == Solve_Succeeded) {
735 CONSOLE_DEBUG("Solution of the primal variables, x");
736 for (i=0; i<sys->n; i++) {
737 CONSOLE_DEBUG(" x[%d] = %e\n", i, x[i]);
738 }
739
740 CONSOLE_DEBUG("Solution of the bound multipliers, z_L and z_U");
741 for (i=0; i<sys->n; i++) {
742 CONSOLE_DEBUG(" z_L[%d] = %e", i, mult_x_L[i]);
743 }
744 for (i=0; i<sys->n; i++) {
745 CONSOLE_DEBUG(" z_U[%d] = %e", i, mult_x_U[i]);
746 }
747
748 CONSOLE_DEBUG("Objective value");
749 CONSOLE_DEBUG(" f(x*) = %e", sys->obj_val);
750
751 ret = 0; /* success */
752 }else{
753 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed solve, unknown status");
754 ret = 1; /* failure */
755 }
756
757 /* free allocated memory */
758 FreeIpoptProblem(sys->nlp);
759 ASC_FREE(x);
760 ASC_FREE(mult_x_L);
761 ASC_FREE(mult_x_U);
762
763 return ret;
764 }
765
766 /**
767 Prepare sys for entering an iteration, increasing the iteration counts
768 and starting the clock.
769 */
770 static void ipopt_iteration_begins(IpoptSystem *sys){
771 sys->clock = tm_cpu_time();
772 ++(sys->s.block.iteration);
773 ++(sys->s.iteration);
774 }
775
776
777 /*
778 Prepare sys for exiting an iteration, stopping the clock and recording
779 the cpu time.
780 */
781 static void ipopt_iteration_ends(IpoptSystem *sys){
782 double cpu_elapsed; /* elapsed this iteration */
783
784 cpu_elapsed = (double)(tm_cpu_time() - sys->clock);
785 sys->s.block.cpu_elapsed += cpu_elapsed;
786 sys->s.cpu_elapsed += cpu_elapsed;
787 }
788
789
790
791 static int ipopt_iterate(slv_system_t server, SlvClientToken asys){
792 UNUSED_PARAMETER(server);
793 CONSOLE_DEBUG("ITERATING...");
794 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
795 return 1;
796 }
797
798 static int ipopt_resolve(slv_system_t server, SlvClientToken asys){
799 UNUSED_PARAMETER(server);
800
801 /* if implementing this, use the 'warm start' thing in IPOPT */
802
803 /* provide initial values of the 'multipliers' */
804
805 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
806 return 1;
807 }
808
809 static const SlvFunctionsT ipopt_internals = {
810 67
811 ,"IPOPT"
812 ,ipopt_create
813 ,ipopt_destroy
814 ,ipopt_eligible_solver
815 ,ipopt_get_default_parameters
816 ,ipopt_get_parameters
817 ,ipopt_set_parameters
818 ,ipopt_get_status
819 ,ipopt_solve
820 ,ipopt_presolve
821 ,ipopt_iterate
822 ,ipopt_resolve
823 ,NULL
824 ,NULL
825 ,NULL
826 };
827
828 int ipopt_register(void){
829 return solver_register(&ipopt_internals);
830 }

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