/[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 1636 - (show annotations) (download) (as text)
Wed Oct 24 13:04:05 2007 UTC (12 years, 3 months ago) by jpye
File MIME type: text/x-csrc
File size: 23773 byte(s)
IPOPT gets installed if built.
added relman_hessian_count placeholder.
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 CONSOLE_DEBUG("Analysing of Hessian matrix sparsity structure not implemented");
488
489 /*
490 for(i=0; i<nvars; ++i){
491 for(j=i; j<nvars; ++j){
492 if(r2(sys->obj, sys->vlist[i], sys->vlist[j]) != 0){
493 iRow[nele_hess] = i;
494 jCol[nele_hess] = j;
495 nele_hess++
496 }
497 }
498 }
499 */
500
501 }else{
502 asc_assert(jCol==NULL);
503 asc_assert(x!=NULL); asc_assert(lambda!=NULL); asc_assert(values!=NULL);
504
505 /* evaluate the Hessian matrix */
506 CONSOLE_DEBUG("Evaluation of Hessian matrix not implemented");
507 }
508
509 return 0; /* fail: not yet implemented */
510 }
511
512 /*------------------------------------------------------------------------------
513 SOLVE ROUTINES
514 */
515
516 static int ipopt_presolve(slv_system_t server, SlvClientToken asys){
517 IpoptSystem *sys;
518 int max, i;
519
520 CONSOLE_DEBUG("PRESOLVE");
521
522 sys = SYS(asys);
523 ipopt_iteration_begins(sys);
524 //check_system(sys);
525
526 asc_assert(sys->vlist && sys->rlist);
527
528 /** @TODO work out if matrix creation is not again needed */
529
530 /** @TODO slv_sort_rels_and_vars(server,&(sys->m),&(sys->n)); */
531
532 /* set all relations as being 'unsatisfied' to start with... */
533 for(i=0; i < sys->rtot; ++i){
534 rel_set_satisfied(sys->rlist[i],FALSE);
535 }
536
537 sys->obj = slv_get_obj_relation(server); /*may have changed objective*/
538 if(!sys->obj){
539 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No objective function was specified");
540 return -3;
541 }
542
543 CONSOLE_DEBUG("got objective rel %p",sys->obj);
544
545 /* calculate nnz for hessian matrix @TODO FIXME */
546
547 sys->nnzH = relman_hessian_count(sys->rlist, sys->rtot, &(sys->vfilt), &(sys->rfilt), &max);
548
549 /* need to provide sparsity structure for hessian matrix? */
550
551 #if 0
552 /** @SEE http://www.coin-or.org/Ipopt/documentation/node37.html */
553 ipopt_eval_h(number_of_variables, NULL/*x at which to evaluate*/, TRUE /* new x */
554 , 1.0/*obj_factor*/, number_of_constraints, lambda/* values of the constraint multipliers */
555 , TRUE /* new lambda */, 0 /* number of nonzero elements in the Hessian */, Index* iRow
556 , Index* jCol, Number* values
557 , void *user_data
558 );
559 #endif
560
561
562 max = relman_obj_direction(sys->obj);
563 if(max==-1){
564 CONSOLE_DEBUG("this is a MINIMIZE problem");
565 }else{
566 CONSOLE_DEBUG("this is a MAXIMIZE problem");
567 }
568
569 CONSOLE_DEBUG("got %d relations and %d vars in system", sys->rtot, sys->vtot);
570 /* calculate number of non-zeros in the Jacobian matrix for the constraint equations */
571
572 sys->nnzJ = relman_jacobian_count(sys->rlist, sys->rtot, &(sys->vfilt), &(sys->rfilt), &max);
573
574 CONSOLE_DEBUG("got %d non-zeros in constraint Jacobian", sys->nnzJ);
575
576 /* need to provide sparsity structure for jacobian? */
577
578
579
580 #if 0
581 if(sys->presolved > 0) { /* system has been presolved before */
582 if(!conopt_dof_changed(sys) /*no changes in fixed or included flags*/
583 && sys->p.partition == sys->J.old_partition
584 && sys->obj == sys->old_obj
585 ){
586 matrix_creation_needed = 0;
587 CONOPT_CONSOLE_DEBUG("YOU JUST AVOIDED MATRIX DESTRUCTION/CREATION");
588 }
589 }
590 #endif
591
592 #if 0
593 // check all this...
594
595 sys->presolved = 1; /* full presolve recognized here */
596 sys->resolve = 0; /* initialize resolve flag */
597
598 sys->J.old_partition = sys->p.partition;
599 sys->old_obj = sys->obj;
600
601 slv_sort_rels_and_vars(server,&(sys->con.m),&(sys->con.n));
602 CONOPT_CONSOLE_DEBUG("FOUND %d CONSTRAINTS AND %d VARS",sys->con.m,sys->con.n);
603 if (sys->obj != NULL) {
604 CONOPT_CONSOLE_DEBUG("ADDING OBJECT AS A ROW");
605 sys->con.m++; /* treat objective as a row */
606 }
607
608 cntvect = ASC_NEW_ARRAY(int,COIDEF_Size());
609 COIDEF_Ini(cntvect);
610 sys->con.cntvect = cntvect;
611 CONOPT_CONSOLE_DEBUG("NUMBER OF CONSTRAINTS = %d",sys->con.m);
612 COIDEF_NumVar(cntvect, &(sys->con.n));
613 COIDEF_NumCon(cntvect, &(sys->con.m));
614 sys->con.nz = num_jacobian_nonzeros(sys, &(sys->con.maxrow));
615 COIDEF_NumNZ(cntvect, &(sys->con.nz));
616 COIDEF_NumNlNz(cntvect, &(sys->con.nz));
617
618 sys->con.base = 0;
619 COIDEF_Base(cntvect,&(sys->con.base));
620 COIDEF_ErrLim(cntvect, &(DOMLIM));
621 COIDEF_ItLim(cntvect, &(ITER_LIMIT));
622
623 if(sys->obj!=NULL){
624 sys->con.optdir = relman_obj_direction(sys->obj);
625 sys->con.objcon = sys->con.m - 1; /* objective will be last row */
626 CONOPT_CONSOLE_DEBUG("SETTING OBJECTIVE CONSTRAINT TO BE %d",sys->con.objcon);
627 }else{
628 sys->con.optdir = 0;
629 sys->con.objcon = 0;
630 }
631 COIDEF_OptDir(cntvect, &(sys->con.optdir));
632 COIDEF_ObjCon(cntvect, &(sys->con.objcon));
633
634 temp = 0;
635 COIDEF_StdOut(cntvect, &temp);
636
637 int debugfv = 1;
638 COIDEF_DebugFV(cntvect, &debugfv);
639
640 destroy_vectors(sys);
641 destroy_matrices(sys);
642 create_matrices(server,sys);
643 create_vectors(sys);
644
645 sys->s.block.current_reordered_block = -2;
646 }
647
648 //...
649
650 if( matrix_creation_needed ) {
651 destroy_array(sys->s.cost);
652 sys->s.cost = create_zero_array(sys->s.costsize,struct slv_block_cost);
653 for( ind = 0; ind < sys->s.costsize; ++ind ) {
654 sys->s.cost[ind].reorder_method = -1;
655 }
656 } else {
657 reset_cost(sys->s.cost,sys->s.costsize);
658 }
659
660 #endif
661
662 /* Reset status */
663 sys->s.iteration = 0;
664 sys->s.cpu_elapsed = 0.0;
665 sys->s.converged = sys->s.diverged = sys->s.inconsistent = FALSE;
666 sys->s.block.previous_total_size = 0;
667 sys->s.costsize = 1+sys->s.block.number_of;
668
669
670 /* set to go to first unconverged block */
671 sys->s.block.current_block = -1;
672 sys->s.block.current_size = 0;
673 sys->s.calc_ok = TRUE;
674 sys->s.block.iteration = 0;
675 sys->obj_val = MAXDOUBLE/2000.0;
676
677 ipopt_iteration_ends(sys);
678
679 CONSOLE_DEBUG("Reset status");
680
681 /* sys->s.cost[sys->s.block.number_of].time=sys->s.cpu_elapsed; */
682
683 ERROR_REPORTER_HERE(ASC_PROG_ERR,"presolve completed");
684 return 0;
685 }
686
687
688 static int ipopt_solve(slv_system_t server, SlvClientToken asys){
689 IpoptSystem *sys;
690 UNUSED_PARAMETER(server);
691 enum ApplicationReturnStatus status;
692 int ret, i, j;
693 struct var_variable *var;
694
695 sys = SYS(asys);
696
697 double *x, *x_L, *x_U, *g_L, *g_U, *mult_x_L, *mult_x_U;
698
699 CONSOLE_DEBUG("SOLVING...");
700
701 /* set the number of variables and allocate space for the bounds */
702 x_L = ASC_NEW_ARRAY(Number,sys->n);
703 x_U = ASC_NEW_ARRAY(Number,sys->n);
704
705 /* @TODO set the values for the variable bounds */
706 int jj = 0;
707 for(j = 0; j < sys->vtot; j++){
708 var = sys->vlist[j];
709 if(var_apply_filter(var,&(sys->vfilt))){
710 x_L[jj] = var_lower_bound(var);
711 x_U[jj] = var_upper_bound(var);
712 jj++;
713 }
714 }
715
716 /** @TODO set bounds on the constraints? */
717 /* need to identify equations that share the same non-constant parts? */
718 /* then find the constant parts and make then g_L or g_U accordingly */
719 /* what to do about other bounds? */
720
721 /* create the IpoptProblem */
722 sys->nlp = CreateIpoptProblem(sys->n, x_L, x_U, sys->m, g_L, g_U, sys->nnzJ, sys->nnzH, 0/*index style=C*/,
723 &ipopt_eval_f, &ipopt_eval_g, &ipopt_eval_grad_f,
724 &ipopt_eval_jac_g, &ipopt_eval_h
725 );
726
727 /* We can free the memory now - the values for the bounds have been
728 copied internally in CreateIpoptProblem */
729 ASC_FREE(x_L);
730 ASC_FREE(x_U);
731 ASC_FREE(g_L);
732 ASC_FREE(g_U);
733
734 /* set some options */
735 AddIpoptNumOption(sys->nlp, "tol", SLV_PARAM_BOOL(&(sys->p),IPOPT_PARAM_TOL));
736 AddIpoptStrOption(sys->nlp, "mu_strategy", SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_MU_STRATEGY));
737
738 /* initial values */
739 x = ASC_NEW_ARRAY(Number, sys->n);
740
741 /** @TODO get values of 'x' from the model */
742
743 /* allocate space to store the bound multipliers at the solution */
744 mult_x_L = ASC_NEW_ARRAY(Number, sys->n);
745 mult_x_U = ASC_NEW_ARRAY(Number, sys->n);
746
747 /* solve the problem */
748 status = IpoptSolve(sys->nlp, x, NULL, &sys->obj_val, NULL, mult_x_L, mult_x_U, NULL);
749
750 /** @TODO update the sys->s.xxxxx flags based on value of 'status' */
751
752 if (status == Solve_Succeeded) {
753 CONSOLE_DEBUG("Solution of the primal variables, x");
754 for (i=0; i<sys->n; i++) {
755 CONSOLE_DEBUG(" x[%d] = %e\n", i, x[i]);
756 }
757
758 CONSOLE_DEBUG("Solution of the bound multipliers, z_L and z_U");
759 for (i=0; i<sys->n; i++) {
760 CONSOLE_DEBUG(" z_L[%d] = %e", i, mult_x_L[i]);
761 }
762 for (i=0; i<sys->n; i++) {
763 CONSOLE_DEBUG(" z_U[%d] = %e", i, mult_x_U[i]);
764 }
765
766 CONSOLE_DEBUG("Objective value");
767 CONSOLE_DEBUG(" f(x*) = %e", sys->obj_val);
768
769 ret = 0; /* success */
770 }else{
771 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed solve, unknown status");
772 ret = 1; /* failure */
773 }
774
775 /* free allocated memory */
776 FreeIpoptProblem(sys->nlp);
777 ASC_FREE(x);
778 ASC_FREE(mult_x_L);
779 ASC_FREE(mult_x_U);
780
781 return ret;
782 }
783
784 /**
785 Prepare sys for entering an iteration, increasing the iteration counts
786 and starting the clock.
787 */
788 static void ipopt_iteration_begins(IpoptSystem *sys){
789 sys->clock = tm_cpu_time();
790 ++(sys->s.block.iteration);
791 ++(sys->s.iteration);
792 }
793
794
795 /*
796 Prepare sys for exiting an iteration, stopping the clock and recording
797 the cpu time.
798 */
799 static void ipopt_iteration_ends(IpoptSystem *sys){
800 double cpu_elapsed; /* elapsed this iteration */
801
802 cpu_elapsed = (double)(tm_cpu_time() - sys->clock);
803 sys->s.block.cpu_elapsed += cpu_elapsed;
804 sys->s.cpu_elapsed += cpu_elapsed;
805 }
806
807
808
809 static int ipopt_iterate(slv_system_t server, SlvClientToken asys){
810 UNUSED_PARAMETER(server);
811 CONSOLE_DEBUG("ITERATING...");
812 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
813 return 1;
814 }
815
816 static int ipopt_resolve(slv_system_t server, SlvClientToken asys){
817 UNUSED_PARAMETER(server);
818
819 /* if implementing this, use the 'warm start' thing in IPOPT */
820
821 /* provide initial values of the 'multipliers' */
822
823 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
824 return 1;
825 }
826
827 static const SlvFunctionsT ipopt_internals = {
828 67
829 ,"IPOPT"
830 ,ipopt_create
831 ,ipopt_destroy
832 ,ipopt_eligible_solver
833 ,ipopt_get_default_parameters
834 ,ipopt_get_parameters
835 ,ipopt_set_parameters
836 ,ipopt_get_status
837 ,ipopt_solve
838 ,ipopt_presolve
839 ,ipopt_iterate
840 ,ipopt_resolve
841 ,NULL
842 ,NULL
843 ,NULL
844 };
845
846 int ipopt_register(void){
847 return solver_register(&ipopt_internals);
848 }

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