/[ascend]/trunk/solvers/ipopt/asc_ipopt.c
ViewVC logotype

Diff of /trunk/solvers/ipopt/asc_ipopt.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1548 by jpye, Mon Jul 23 06:25:49 2007 UTC revision 1549 by jpye, Mon Jul 23 14:30:35 2007 UTC
# Line 64  ASC_DLLSPEC SolverRegisterFn ipopt_regis Line 64  ASC_DLLSPEC SolverRegisterFn ipopt_regis
64  enum{  enum{
65      IPOPT_PARAM_TOL      IPOPT_PARAM_TOL
66      ,IPOPT_PARAM_MAX_ITER      ,IPOPT_PARAM_MAX_ITER
67      ,IPOPT_PARAM_SAFECALC      ,IPOPT_PARAM_SAFEEVAL
68      ,IPOPT_PARAM_MU_STRATEGY      ,IPOPT_PARAM_MU_STRATEGY
69      ,IPOPT_PARAMS      ,IPOPT_PARAMS
70  };  };
# Line 81  struct IpoptSystemStruct{ Line 81  struct IpoptSystemStruct{
81      struct rel_relation         *old_obj;/* Objective function: NULL = none */      struct rel_relation         *old_obj;/* Objective function: NULL = none */
82      struct var_variable         **vlist; /* Variable list (NULL terminated) */      struct var_variable         **vlist; /* Variable list (NULL terminated) */
83      struct rel_relation         **rlist; /* Relation list (NULL terminated) */      struct rel_relation         **rlist; /* Relation list (NULL terminated) */
84        var_filter_t vfilt;
85    
86      /*      /*
87          Solver information          Solver information
# Line 99  struct IpoptSystemStruct{ Line 100  struct IpoptSystemStruct{
100      double                 clock;        /* CPU time */      double                 clock;        /* CPU time */
101    
102      int32 calc_ok;      int32 calc_ok;
103        double obj_val;
104    
105      void *parm_array[IPOPT_PARAMS];      void *parm_array[IPOPT_PARAMS];
106      struct slv_parameter pa[IPOPT_PARAMS];      struct slv_parameter pa[IPOPT_PARAMS];
# Line 110  struct IpoptSystemStruct{ Line 112  struct IpoptSystemStruct{
112      Index m;                          /* number of constraints */      Index m;                          /* number of constraints */
113    
114      int nnzJ; /* number of non zeros in the jacobian of the constraints */      int nnzJ; /* number of non zeros in the jacobian of the constraints */
115        int nnzH; /* number of non-zeros in the hessian of the objective */
116    
117      Number* x_L;                  /* lower bounds on x */      Number* x_L;                  /* lower bounds on x */
118      Number* x_U;                  /* upper bounds on x */      Number* x_U;                  /* upper bounds on x */
# Line 129  typedef struct IpoptSystemStruct IpoptSy Line 132  typedef struct IpoptSystemStruct IpoptSy
132    
133  static int ipopt_get_default_parameters(slv_system_t server, SlvClientToken asys, slv_parameters_t *parameters);  static int ipopt_get_default_parameters(slv_system_t server, SlvClientToken asys, slv_parameters_t *parameters);
134    
135    static void ipopt_iteration_begins(IpoptSystem *sys);
136    static void ipopt_iteration_ends(IpoptSystem *sys);
137    
138  /*------------------------------------------------------------------------------  /*------------------------------------------------------------------------------
139    SYSTEM SETUP/DESTROY AND SOLVER ELIGIBILITY    SYSTEM SETUP/DESTROY AND SOLVER ELIGIBILITY
140  */  */
# Line 157  static SlvClientToken ipopt_create(slv_s Line 163  static SlvClientToken ipopt_create(slv_s
163      sys->s.calc_ok = TRUE;      sys->s.calc_ok = TRUE;
164      sys->s.costsize = 0;      sys->s.costsize = 0;
165      sys->s.cost = NULL; /*redundant, but sanity preserving */      sys->s.cost = NULL; /*redundant, but sanity preserving */
166        sys->s.block.number_of = 1;
167        sys->s.block.current_block = 0;
168        sys->s.block.current_reordered_block = 0;
169        sys->s.block.current_size = 0;
170        sys->s.block.previous_total_size = 0;
171        sys->s.block.iteration = 0;
172        sys->s.block.funcs = 0;
173        sys->s.block.jacs = 0;
174        sys->s.block.cpu_elapsed = 0;
175        sys->s.block.functime = 0;
176        sys->s.block.jactime = 0;
177        sys->s.block.residual = 0;
178    
179        /** @TODO set up sys->vfilt */
180    
181      sys->vlist = slv_get_solvers_var_list(server);      sys->vlist = slv_get_solvers_var_list(server);
182      sys->rlist = slv_get_solvers_rel_list(server);      sys->rlist = slv_get_solvers_rel_list(server);
# Line 183  static SlvClientToken ipopt_create(slv_s Line 203  static SlvClientToken ipopt_create(slv_s
203    
204  static int32 ipopt_destroy(slv_system_t server, SlvClientToken asys){  static int32 ipopt_destroy(slv_system_t server, SlvClientToken asys){
205      UNUSED_PARAMETER(server);      UNUSED_PARAMETER(server);
206      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"ipopt_destroy not implemented");
207      return 1;      return 1;
208  }    }  
209    
210  static int32 ipopt_eligible_solver(slv_system_t server){  static int32 ipopt_eligible_solver(slv_system_t server){
211      UNUSED_PARAMETER(server);      UNUSED_PARAMETER(server);
212      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");      
213        /// TODO check that there is a MAXIMIZE or MINIMIZE statement
214        /// TODO check that there are no discrete-valued variables
215        /// TODO check that there are no WHENs or CONDITIONALs
216        /// TODO check anything else?
217    
218        ERROR_REPORTER_HERE(ASC_PROG_WARNING,"ipopt_eligible_solver not implemented");
219      return 1;      return 1;
220  }  }
221    
# Line 246  int32 ipopt_get_default_parameters(slv_s Line 272  int32 ipopt_get_default_parameters(slv_s
272          ,(SlvParameterInitChar){{"mu_strategy"          ,(SlvParameterInitChar){{"mu_strategy"
273              ,"Update strategy for barrier parameter",6              ,"Update strategy for barrier parameter",6
274              ,"Determines which barrier parameter update strategy is to be used."              ,"Determines which barrier parameter update strategy is to be used."
275              " 'MONOTONE' is the monotone (Fiacco-McCormick) strategy;"              " 'monotone' is the monotone (Fiacco-McCormick) strategy;"
276              " 'ADAPTIVE' is the adaptive update strategy."              " 'adaptive' is the adaptive update strategy."
277          }, "MONOTONE"}, (char *[]){          }, "monotone"}, (char *[]){
278              "MONOTONE","ADAPTIVE",NULL              "monotone","adaptive",NULL
279          }          }
280      );      );
281    
282        slv_param_bool(parameters,IPOPT_PARAM_SAFEEVAL
283            ,(SlvParameterInitBool){{"safeeval"
284                ,"Use safe evaluation?",1
285                ,"Use 'safe' function evaluation routines (TRUE) or allow ASCEND to "
286                "throw SIGFPE errors which will then halt integration (FALSE)."
287            }, FALSE}
288        );
289    
290      asc_assert(parameters->num_parms==IPOPT_PARAMS);      asc_assert(parameters->num_parms==IPOPT_PARAMS);
291    
292      return 1;      return 1;
# Line 282  static void ipopt_set_parameters(slv_sys Line 316  static void ipopt_set_parameters(slv_sys
316  }  }
317    
318  /*------------------------------------------------------------------------------  /*------------------------------------------------------------------------------
319    EVALUATION FUNCTIONS    EVALUATION AND RESULT HOOK FUNCTIONS
320  */  */
321    
322  /**  /**
# Line 317  Bool ipopt_eval_f(Index n, Number *x, Bo Line 351  Bool ipopt_eval_f(Index n, Number *x, Bo
351    
352      sys->calc_ok = TRUE;      sys->calc_ok = TRUE;
353    
354      *obj_value = relman_eval(sys->obj,&(sys->calc_ok),SLV_PARAM_BOOL(&(sys->p),IPOPT_PARAM_SAFECALC));      *obj_value = relman_eval(sys->obj,&(sys->calc_ok),SLV_PARAM_BOOL(&(sys->p),IPOPT_PARAM_SAFEEVAL));
355    
356      return sys->calc_ok;      return sys->calc_ok;
357  }  }
# Line 357  Bool ipopt_eval_grad_f(Index n, Number* Line 391  Bool ipopt_eval_grad_f(Index n, Number*
391    
392      relman_diff2(      relman_diff2(
393          sys->obj,&vfilter,derivatives,variables          sys->obj,&vfilter,derivatives,variables
394          , &count,SLV_PARAM_BOOL(&(sys->p),IPOPT_PARAM_SAFECALC)          , &count,SLV_PARAM_BOOL(&(sys->p),IPOPT_PARAM_SAFEEVAL)
395      );      );
396    
397      for(j=0; j<len; ++j){      for(j=0; j<len; ++j){
# Line 396  Bool ipopt_eval_jac_g(Index n, Number* x Line 430  Bool ipopt_eval_jac_g(Index n, Number* x
430  ){  ){
431      IpoptSystem *sys;      IpoptSystem *sys;
432      sys = SYS(user_data);      sys = SYS(user_data);
433      int i,j,res;      int i,res;
434    
435      asc_assert(n==sys->n);      asc_assert(n==sys->n);
436      asc_assert(nele_jac==sys->nnzJ);      asc_assert(nele_jac==sys->nnzJ);
# Line 420  Bool ipopt_eval_h(Index n, Number* x, Bo Line 454  Bool ipopt_eval_h(Index n, Number* x, Bo
454          , Number obj_factor, Index m, Number* lambda          , Number obj_factor, Index m, Number* lambda
455          , Bool new_lambda, Index nele_hess, Index* iRow          , Bool new_lambda, Index nele_hess, Index* iRow
456          , Index* jCol, Number* values          , Index* jCol, Number* values
457            , void *user_data
458  ){  ){
459      /* not sure about this one yet: do all the 2nd deriv things work for this? */      /* not sure about this one yet: do all the 2nd deriv things work for this? */
460      return 0; /* fail: not yet implemented */      return 0; /* fail: not yet implemented */
# Line 432  Bool ipopt_eval_h(Index n, Number* x, Bo Line 467  Bool ipopt_eval_h(Index n, Number* x, Bo
467  static int ipopt_solve(slv_system_t server, SlvClientToken asys){  static int ipopt_solve(slv_system_t server, SlvClientToken asys){
468      IpoptSystem *sys;      IpoptSystem *sys;
469      UNUSED_PARAMETER(server);      UNUSED_PARAMETER(server);
470      int i;      enum ApplicationReturnStatus status;
471        int ret, i, j;
472        struct var_variable *var;
473    
474      sys = SYS(asys);      sys = SYS(asys);
475    
476      /* set the number of variables and allocate space for the bounds */      double *x, *x_L, *x_U, *g_L, *g_U, *mult_x_L, *mult_x_U;
     sys->x_L = ASC_NEW_ARRAY(Number,sys->n);  
     sys->x_U = ASC_NEW_ARRAY(Number,sys->n);  
477    
478      /* set the values for the variable bounds */      /* set the number of variables and allocate space for the bounds */
479        x_L = ASC_NEW_ARRAY(Number,sys->n);
480        x_U = ASC_NEW_ARRAY(Number,sys->n);
481    
482  #if 0      /* @TODO set the values for the variable bounds */
483  // sample code for the C interface      int jj = 0;
484        for(j = 0; j < sys->vtot; j++){
485            var = sys->vlist[j];
486            if(var_apply_filter(var,&(sys->vfilt))){
487                x_L[jj] = var_lower_bound(var);
488                x_U[jj] = var_upper_bound(var);
489                jj++;
490            }
491        }
492    
493  int main(){      /** @TODO set bounds on the constraints? */
494    
495    /* set the number of constraints and allocate space for the bounds */      /* create the IpoptProblem */
496    m=2;      sys->nlp = CreateIpoptProblem(sys->n, x_L, x_U, sys->m, g_L, g_U, sys->nnzJ, sys->nnzH, 0/*index style=C*/,
497    g_L = (Number*)malloc(sizeof(Number)*m);          &ipopt_eval_f, &ipopt_eval_g, &ipopt_eval_grad_f,
498    g_U = (Number*)malloc(sizeof(Number)*m);          &ipopt_eval_jac_g, &ipopt_eval_h
499    /* set the values of the constraint bounds */      );
   g_L[0] = 25; g_U[0] = 2e19;  
   g_L[1] = 40; g_U[1] = 40;  
   
   /* create the IpoptProblem */  
   nlp = CreateIpoptProblem(n, x_L, x_U, m, g_L, g_U, 8, 10, 0,  
                &eval_f, &eval_g, &eval_grad_f,  
                &eval_jac_g, &eval_h);  
500        
501    /* We can free the memory now - the values for the bounds have been      /* We can free the memory now - the values for the bounds have been
502       copied internally in CreateIpoptProblem */      copied internally in CreateIpoptProblem */
503    free(x_L);      ASC_FREE(x_L);
504    free(x_U);      ASC_FREE(x_U);
505    free(g_L);      ASC_FREE(g_L);
506    free(g_U);      ASC_FREE(g_U);
507    
508    /* set some options */      /* set some options */
509    AddIpoptNumOption(nlp, "tol", 1e-9);      AddIpoptNumOption(sys->nlp, "tol", SLV_PARAM_BOOL(&(sys->p),IPOPT_PARAM_TOL));
510    AddIpoptStrOption(nlp, "mu_strategy", "adaptive");      AddIpoptStrOption(sys->nlp, "mu_strategy", SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_MU_STRATEGY));
511    
512    /* allocate space for the initial point and set the values */      /* initial values */
513    x = (Number*)malloc(sizeof(Number)*n);      x = ASC_NEW_ARRAY(Number, sys->n);
514    x[0] = 1.0;  
515    x[1] = 5.0;      /** @TODO get values of 'x' from the model */
516    x[2] = 5.0;  
517    x[3] = 1.0;      /* allocate space to store the bound multipliers at the solution */
518        mult_x_L = ASC_NEW_ARRAY(Number, sys->n);
519    /* allocate space to store the bound multipliers at the solution */      mult_x_U = ASC_NEW_ARRAY(Number, sys->n);
520    mult_x_L = (Number*)malloc(sizeof(Number)*n);  
521    mult_x_U = (Number*)malloc(sizeof(Number)*n);      /* solve the problem */
522        status = IpoptSolve(sys->nlp, x, NULL, &sys->obj_val, NULL, mult_x_L, mult_x_U, NULL);
523    /* solve the problem */  
524    status = IpoptSolve(nlp, x, NULL, &obj, NULL, mult_x_L, mult_x_U, NULL);      /** @TODO update the sys->s.xxxxx flags based on value of 'status' */
525    
526    if (status == Solve_Succeeded) {      if (status == Solve_Succeeded) {
527      printf("\n\nSolution of the primal variables, x\n");          CONSOLE_DEBUG("Solution of the primal variables, x");
528      for (i=0; i<n; i++) {          for (i=0; i<sys->n; i++) {
529        printf("x[%d] = %e\n", i, x[i]);              CONSOLE_DEBUG("   x[%d] = %e\n", i, x[i]);
530      }          }
   
     printf("\n\nSolution of the bound multipliers, z_L and z_U\n");  
     for (i=0; i<n; i++) {  
       printf("z_L[%d] = %e\n", i, mult_x_L[i]);  
     }  
     for (i=0; i<n; i++) {  
       printf("z_U[%d] = %e\n", i, mult_x_U[i]);  
     }  
   
     printf("\n\nObjective value\n");  
     printf("f(x*) = %e\n", obj);  
   }  
   
   /* free allocated memory */  
   FreeIpoptProblem(nlp);  
   free(x);  
   free(mult_x_L);  
   free(mult_x_U);  
531    
532    return 0;          CONSOLE_DEBUG("Solution of the bound multipliers, z_L and z_U");
533  }          for (i=0; i<sys->n; i++) {
534  #endif              CONSOLE_DEBUG("   z_L[%d] = %e", i, mult_x_L[i]);
535            }
536            for (i=0; i<sys->n; i++) {
537                CONSOLE_DEBUG("    z_U[%d] = %e", i, mult_x_U[i]);
538            }
539    
540            CONSOLE_DEBUG("Objective value");
541            CONSOLE_DEBUG("    f(x*) = %e", sys->obj_val);
542    
543      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");          ret = 0; /* success */
544      return 1;      }else{
545            ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed solve, unknown status");
546            ret = 1; /* failure */
547        }
548    
549        /* free allocated memory */
550        FreeIpoptProblem(sys->nlp);
551        ASC_FREE(x);
552        ASC_FREE(mult_x_L);
553        ASC_FREE(mult_x_U);
554    
555        return ret;
556  }  }
557    
558  static int ipopt_presolve(slv_system_t server, SlvClientToken asys){  static int ipopt_presolve(slv_system_t server, SlvClientToken asys){
559      IpoptSystem *sys;      IpoptSystem *sys;
     struct var_variable **vp;  
     struct rel_relation **rp;  
     int cap, ind;  
     int matrix_creation_needed = 1;  
     int *cntvect, temp;  
560    
561      CONSOLE_DEBUG("PRESOLVE");      CONSOLE_DEBUG("PRESOLVE");
562    
563      sys = SYS(asys);      sys = SYS(asys);
564      //iteration_begins(sys);      ipopt_iteration_begins(sys);
565      //check_system(sys);      //check_system(sys);
566    
567      asc_assert(sys->vlist && sys->rlist);      asc_assert(sys->vlist && sys->rlist);
# Line 579  static int ipopt_presolve(slv_system_t s Line 611  static int ipopt_presolve(slv_system_t s
611          sys->presolved = 1; /* full presolve recognized here */          sys->presolved = 1; /* full presolve recognized here */
612          sys->resolve = 0;   /* initialize resolve flag */          sys->resolve = 0;   /* initialize resolve flag */
613    
614          /* provide nnz for jacobian matrix of constraints */          /* @TODO calculate nnz for jacobian matrix of constraints */
615    
616          /* provide sparsity structure for jacobian */          /* provide sparsity structure for jacobian */
617    
618          /* provide nnz for hessian matrix */          /** @TODO calculate nnz for hessian matrix */
619    
620          /* provide sparsity structure for hessian matrix */          /* provide sparsity structure for hessian matrix */
621    
# Line 646  static int ipopt_presolve(slv_system_t s Line 678  static int ipopt_presolve(slv_system_t s
678          sys->s.block.current_reordered_block = -2;          sys->s.block.current_reordered_block = -2;
679      }      }
680    
681      /* Reset status */      //...
     sys->con.optimized = 0;  
     sys->s.iteration = 0;  
     sys->s.cpu_elapsed = 0.0;  
     sys->s.converged = sys->s.diverged = sys->s.inconsistent = FALSE;  
     sys->s.block.previous_total_size = 0;  
     sys->s.costsize = 1+sys->s.block.number_of;  
682    
683      if( matrix_creation_needed ) {      if( matrix_creation_needed ) {
684          destroy_array(sys->s.cost);          destroy_array(sys->s.cost);
# Line 664  static int ipopt_presolve(slv_system_t s Line 690  static int ipopt_presolve(slv_system_t s
690          reset_cost(sys->s.cost,sys->s.costsize);          reset_cost(sys->s.cost,sys->s.costsize);
691      }      }
692    
693    #endif
694    
695        /* Reset status */
696        sys->s.iteration = 0;
697        sys->s.cpu_elapsed = 0.0;
698        sys->s.converged = sys->s.diverged = sys->s.inconsistent = FALSE;
699        sys->s.block.previous_total_size = 0;
700        sys->s.costsize = 1+sys->s.block.number_of;
701    
702      /* set to go to first unconverged block */      /* set to go to first unconverged block */
703      sys->s.block.current_block = -1;      sys->s.block.current_block = -1;
704      sys->s.block.current_size = 0;      sys->s.block.current_size = 0;
705      sys->s.calc_ok = TRUE;      sys->s.calc_ok = TRUE;
706      sys->s.block.iteration = 0;      sys->s.block.iteration = 0;
707      sys->objective =  MAXDOUBLE/2000.0;      sys->obj_val =  MAXDOUBLE/2000.0;
708    
709      update_status(sys);      ipopt_iteration_ends(sys);
     iteration_ends(sys);  
710      sys->s.cost[sys->s.block.number_of].time=sys->s.cpu_elapsed;      sys->s.cost[sys->s.block.number_of].time=sys->s.cpu_elapsed;
 #endif  
711    
712        ERROR_REPORTER_HERE(ASC_PROG_ERR,"presolve completed");
713      return 0;      return 0;
714  }  }
715    
716    /**
717        Prepare sys for entering an iteration, increasing the iteration counts
718        and starting the clock.
719    */
720    static void ipopt_iteration_begins(IpoptSystem *sys){
721        sys->clock = tm_cpu_time();
722        ++(sys->s.block.iteration);
723        ++(sys->s.iteration);
724    }
725    
726    
727    /*
728        Prepare sys for exiting an iteration, stopping the clock and recording
729        the cpu time.
730    */
731    static void ipopt_iteration_ends(IpoptSystem *sys){
732        double cpu_elapsed;   /* elapsed this iteration */
733    
734        cpu_elapsed = (double)(tm_cpu_time() - sys->clock);
735        sys->s.block.cpu_elapsed += cpu_elapsed;
736        sys->s.cpu_elapsed += cpu_elapsed;
737    }
738    
739    
740    
741  static int ipopt_iterate(slv_system_t server, SlvClientToken asys){  static int ipopt_iterate(slv_system_t server, SlvClientToken asys){
742      UNUSED_PARAMETER(server);      UNUSED_PARAMETER(server);
743      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");

Legend:
Removed from v.1548  
changed lines
  Added in v.1549

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