/[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 1554 by jpye, Tue Jul 24 08:56:59 2007 UTC revision 1555 by jpye, Tue Jul 24 14:16:45 2007 UTC
# 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    
85      var_filter_t vfilt;      var_filter_t vfilt;
86        rel_filter_t rfilt;
87    
88      /*      /*
89          Solver information          Solver information
# Line 176  static SlvClientToken ipopt_create(slv_s Line 178  static SlvClientToken ipopt_create(slv_s
178      sys->s.block.jactime = 0;      sys->s.block.jactime = 0;
179      sys->s.block.residual = 0;      sys->s.block.residual = 0;
180    
181      /** @TODO set up sys->vfilt */      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);      sys->vlist = slv_get_solvers_var_list(server);
187      sys->rlist = slv_get_solvers_rel_list(server);      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);      sys->obj = slv_get_obj_relation(server);
193      if(sys->vlist == NULL) {      if(sys->vlist == NULL) {
194          ASC_FREE(sys);          ASC_FREE(sys);
# Line 469  Bool ipopt_eval_h(Index n, Number* x, Bo Line 478  Bool ipopt_eval_h(Index n, Number* x, Bo
478          , Index* jCol, Number* values          , Index* jCol, Number* values
479          , void *user_data          , void *user_data
480  ){  ){
481      /* not sure about this one yet: do all the 2nd deriv things work for this? */      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 */      return 0; /* fail: not yet implemented */
508  }  }
509    
# Line 477  Bool ipopt_eval_h(Index n, Number* x, Bo Line 511  Bool ipopt_eval_h(Index n, Number* x, Bo
511    SOLVE ROUTINES    SOLVE ROUTINES
512  */  */
513    
514  static int ipopt_solve(slv_system_t server, SlvClientToken asys){  static int ipopt_presolve(slv_system_t server, SlvClientToken asys){
515      IpoptSystem *sys;      IpoptSystem *sys;
516      UNUSED_PARAMETER(server);      int max, i;
517      enum ApplicationReturnStatus status;  
518      int ret, i, j;      CONSOLE_DEBUG("PRESOLVE");
     struct var_variable *var;  
519    
520      sys = SYS(asys);      sys = SYS(asys);
521        ipopt_iteration_begins(sys);
522        //check_system(sys);
523    
524      double *x, *x_L, *x_U, *g_L, *g_U, *mult_x_L, *mult_x_U;      asc_assert(sys->vlist && sys->rlist);
525    
526      CONSOLE_DEBUG("SOLVING...");      /** @TODO work out if matrix creation is not again needed */
527    
528      /* set the number of variables and allocate space for the bounds */      /** @TODO slv_sort_rels_and_vars(server,&(sys->m),&(sys->n)); */
     x_L = ASC_NEW_ARRAY(Number,sys->n);  
     x_U = ASC_NEW_ARRAY(Number,sys->n);  
529    
530      /* @TODO set the values for the variable bounds */      /* set all relations as being 'unsatisfied' to start with... */
531      int jj = 0;      for(i=0; i < sys->rtot; ++i){
532      for(j = 0; j < sys->vtot; j++){          rel_set_satisfied(sys->rlist[i],FALSE);
         var = sys->vlist[j];  
         if(var_apply_filter(var,&(sys->vfilt))){  
             x_L[jj] = var_lower_bound(var);  
             x_U[jj] = var_upper_bound(var);  
             jj++;  
         }  
533      }      }
534    
535      /** @TODO set bounds on the constraints? */      sys->obj = slv_get_obj_relation(server); /*may have changed objective*/
536        if(!sys->obj){
537      /* create the IpoptProblem */          ERROR_REPORTER_HERE(ASC_PROG_ERR,"No objective function was specified");
538      sys->nlp = CreateIpoptProblem(sys->n, x_L, x_U, sys->m, g_L, g_U, sys->nnzJ, sys->nnzH, 0/*index style=C*/,          return -3;
539          &ipopt_eval_f, &ipopt_eval_g, &ipopt_eval_grad_f,      }
         &ipopt_eval_jac_g, &ipopt_eval_h  
     );  
     
     /* We can free the memory now - the values for the bounds have been  
     copied internally in CreateIpoptProblem */  
     ASC_FREE(x_L);  
     ASC_FREE(x_U);  
     ASC_FREE(g_L);  
     ASC_FREE(g_U);  
   
     /* set some options */  
     AddIpoptNumOption(sys->nlp, "tol", SLV_PARAM_BOOL(&(sys->p),IPOPT_PARAM_TOL));  
     AddIpoptStrOption(sys->nlp, "mu_strategy", SLV_PARAM_CHAR(&(sys->p),IPOPT_PARAM_MU_STRATEGY));  
   
     /* initial values */  
     x = ASC_NEW_ARRAY(Number, sys->n);  
   
     /** @TODO get values of 'x' from the model */  
   
     /* allocate space to store the bound multipliers at the solution */  
     mult_x_L = ASC_NEW_ARRAY(Number, sys->n);  
     mult_x_U = ASC_NEW_ARRAY(Number, sys->n);  
   
     /* solve the problem */  
     status = IpoptSolve(sys->nlp, x, NULL, &sys->obj_val, NULL, mult_x_L, mult_x_U, NULL);  
540    
541      /** @TODO update the sys->s.xxxxx flags based on value of 'status' */      CONSOLE_DEBUG("got objective rel %p",sys->obj);
542    
543      if (status == Solve_Succeeded) {      /** @TODO calculate nnz for hessian matrix */
         CONSOLE_DEBUG("Solution of the primal variables, x");  
         for (i=0; i<sys->n; i++) {  
             CONSOLE_DEBUG("   x[%d] = %e\n", i, x[i]);  
         }  
544    
545          CONSOLE_DEBUG("Solution of the bound multipliers, z_L and z_U");      /* need to provide sparsity structure for hessian matrix? */
         for (i=0; i<sys->n; i++) {  
             CONSOLE_DEBUG("   z_L[%d] = %e", i, mult_x_L[i]);  
         }  
         for (i=0; i<sys->n; i++) {  
             CONSOLE_DEBUG("    z_U[%d] = %e", i, mult_x_U[i]);  
         }  
546    
547          CONSOLE_DEBUG("Objective value");      max = relman_obj_direction(sys->obj);
548          CONSOLE_DEBUG("    f(x*) = %e", sys->obj_val);      if(max==-1){
549            CONSOLE_DEBUG("this is a MINIMIZE problem");
         ret = 0; /* success */  
550      }else{      }else{
551          ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed solve, unknown status");          CONSOLE_DEBUG("this is a MAXIMIZE problem");
         ret = 1; /* failure */  
552      }      }
   
     /* free allocated memory */  
     FreeIpoptProblem(sys->nlp);  
     ASC_FREE(x);  
     ASC_FREE(mult_x_L);  
     ASC_FREE(mult_x_U);  
553    
554      return ret;      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  static int ipopt_presolve(slv_system_t server, SlvClientToken asys){      sys->nnzJ = relman_jacobian_count(sys->rlist, sys->rtot, &(sys->vfilt), &(sys->rfilt), &max);
     IpoptSystem *sys;  
558    
559      CONSOLE_DEBUG("PRESOLVE");      CONSOLE_DEBUG("got %d non-zeros in constraint Jacobian", sys->nnzJ);
560        
561        /* need to provide sparsity structure for jacobian? */
562    
     sys = SYS(asys);  
     ipopt_iteration_begins(sys);  
     //check_system(sys);  
563    
     asc_assert(sys->vlist && sys->rlist);  
   
     sys->obj = slv_get_obj_relation(server); /*may have changed objective*/  
     if(!sys->obj){  
         ERROR_REPORTER_HERE(ASC_PROG_ERR,"No objective function was specified");  
         return -3;  
     }  
   
     CONSOLE_DEBUG("got objective rel %p",sys->obj);  
564    
565  #if 0  #if 0
566      if(sys->presolved > 0) { /* system has been presolved before */      if(sys->presolved > 0) { /* system has been presolved before */
# Line 604  static int ipopt_presolve(slv_system_t s Line 577  static int ipopt_presolve(slv_system_t s
577  #if 0  #if 0
578      // check all this...      // check all this...
579    
     /* set all relations as being 'unsatisfied' to start with... */  
     rp=sys->rlist;  
     for( ind = 0; ind < sys->rtot; ++ind ) {  
         rel_set_satisfied(rp[ind],FALSE);  
     }  
   
     if( matrix_creation_needed ) {  
   
         cap = slv_get_num_solvers_rels(server);  
         sys->cap = slv_get_num_solvers_vars(server);  
         sys->cap = MAX(sys->cap,cap);  
         vp=sys->vlist;  
         for( ind = 0; ind < sys->vtot; ++ind ) {  
             var_set_in_block(vp[ind],FALSE);  
         }  
         rp=sys->rlist;  
         for( ind = 0; ind < sys->rtot; ++ind ) {  
             rel_set_in_block(rp[ind],FALSE);  
             /* rel_set_satisfied(rp[ind],FALSE); */  
         }  
   
580          sys->presolved = 1; /* full presolve recognized here */          sys->presolved = 1; /* full presolve recognized here */
581          sys->resolve = 0;   /* initialize resolve flag */          sys->resolve = 0;   /* initialize resolve flag */
582    
         /* @TODO calculate nnz for jacobian matrix of constraints */  
   
         /* provide sparsity structure for jacobian */  
   
         /** @TODO calculate nnz for hessian matrix */  
   
         /* provide sparsity structure for hessian matrix */  
   
583          sys->J.old_partition = sys->p.partition;          sys->J.old_partition = sys->p.partition;
584          sys->old_obj = sys->obj;          sys->old_obj = sys->obj;
585    
# Line 675  static int ipopt_presolve(slv_system_t s Line 619  static int ipopt_presolve(slv_system_t s
619          temp = 0;          temp = 0;
620          COIDEF_StdOut(cntvect, &temp);          COIDEF_StdOut(cntvect, &temp);
621    
         COIDEF_ReadMatrix(cntvect, &conopt_readmatrix);  
         COIDEF_FDEval(cntvect, &conopt_fdeval);  
         COIDEF_Option(cntvect, &conopt_option);  
         COIDEF_Solution(cntvect, &conopt_solution);  
         COIDEF_Status(cntvect, &conopt_status);  
         COIDEF_Message(cntvect, &asc_conopt_message);  
         COIDEF_ErrMsg(cntvect, &conopt_errmsg);  
         COIDEF_Progress(cntvect, &asc_conopt_progress);  
   
622          int debugfv = 1;          int debugfv = 1;
623          COIDEF_DebugFV(cntvect, &debugfv);          COIDEF_DebugFV(cntvect, &debugfv);
624    
# Line 734  static int ipopt_presolve(slv_system_t s Line 669  static int ipopt_presolve(slv_system_t s
669      return 0;      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      Prepare sys for entering an iteration, increasing the iteration counts
768      and starting the clock.      and starting the clock.

Legend:
Removed from v.1554  
changed lines
  Added in v.1555

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