solvers/ipopt/asc_ipopt.c

/*  ASCEND modelling environment
    Copyright (C) 2007 Carnegie Mellon University

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2, or (at your option)
    any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place - Suite 330,
    Boston, MA 02111-1307, USA.
*//**
    @file
    Connection of the IPOPT optimisation solver into ASCEND.

    THIS IS STILL VERY MUCH UNDER DEVELOPMENT AND INCOMPLETE. I'VE ACTUALLY
    ONLY JUST STARTED WRITING IT by starting with asc_tron.c and modifying.

    The IPOPT solver is documented at http://projects.coin-or.org/Ipopt/
*//*
    ASCEND wrapper for IPOPT originally by John Pye, Jun 2007.
*/

#include <utilities/config.h>

#ifndef WITH_IPOPT
# error "WITH_IPOPT must be defined in order to build this." 
#endif

#include <solver/solver.h>

#include <system/calc.h>
#include <system/relman.h>
#include <system/slv_stdcalls.h>
#include <system/block.h>

#include <utilities/ascConfig.h>
#include <utilities/ascPanic.h>
#include <utilities/ascMalloc.h>
#include <utilities/ascDynaLoad.h>
#include <utilities/mem.h>
#include <utilities/ascEnvVar.h>
#include <general/tm_time.h>
#include <general/env.h>

#include <ipopt/IpStdCInterface.h>

ASC_DLLSPEC SolverRegisterFn ipopt_register;

/*------------------------------------------------------------------------------
  DATA STRUCTURES AND FORWARD DEFS
*/

/**
    Documentation of solver options for IPOPT is at
    http://www.coin-or.org/Ipopt/documentation/node1.html
*/
enum{
    IPOPT_PARAM_TOL
    ,IPOPT_PARAM_MAX_ITER
    ,IPOPT_PARAM_SAFECALC
    ,IPOPT_PARAM_MU_STRATEGY
    ,IPOPT_PARAMS
};

#define SYS(s) ((IpoptSystem *)(s))

struct IpoptSystemStruct{

    /*
        Problem definition
    */
    slv_system_t                slv;     /* slv_system_t back-link */
    struct rel_relation         *obj;    /* Objective function: NULL = none */
    struct rel_relation         *old_obj;/* Objective function: NULL = none */
    struct var_variable         **vlist; /* Variable list (NULL terminated) */
    struct rel_relation         **rlist; /* Relation list (NULL terminated) */

    /*
        Solver information
    */
    int32                  presolved;    /* ? Has the system been presolved */
    int32                  resolve;      /* ? Has the system been resolved */
    slv_parameters_t       p;            /* Parameters */
    slv_status_t           s;            /* Status (as of iteration end) */

    int32                  cap;          /* Order of matrix/vectors */
    int32                  rank;         /* Symbolic rank of problem */
    int32                  vused;        /* Free and incident variables */
    int32                  vtot;         /* length of varlist */
    int32                  rused;        /* Included relations */
    int32                  rtot;         /* length of rellist */
    double                 clock;        /* CPU time */

    void *parm_array[IPOPT_PARAMS];
    struct slv_parameter pa[IPOPT_PARAMS];

    /*
        IPOPT DATA
    */
    Index n;                          /* number of variables */
    Index m;                          /* number of constraints */

    int nnzJ; /* number of non zeros in the jacobian of the constraints */

    Number* x_L;                  /* lower bounds on x */
    Number* x_U;                  /* upper bounds on x */
    Number* g_L;                  /* lower bounds on g */
    Number* g_U;                  /* upper bounds on g */

    IpoptProblem nlp;             /* IpoptProblem */

    enum ApplicationReturnStatus status; /* Solve return code */
    Number* x;                    /* starting point and solution vector */
    Number* mult_x_L;             /* lower bound multipliers at the solution */
    Number* mult_x_U;             /* upper bound multipliers at the solution */
    Number objval;                /* objective value */
    Index i;                      /* generic counter */
};

typedef struct IpoptSystemStruct IpoptSystem;

static int ipopt_get_default_parameters(slv_system_t server, SlvClientToken asys, slv_parameters_t *parameters);

/*------------------------------------------------------------------------------
  SYSTEM SETUP/DESTROY AND SOLVER ELIGIBILITY
*/

static SlvClientToken ipopt_create(slv_system_t server, int32*statusindex){
    IpoptSystem *sys;

    sys = ASC_NEW_CLEAR(IpoptSystem);
    if(sys==NULL){
        *statusindex = 1;
        return sys;
    }

    sys->p.parms = sys->pa;
    sys->p.dynamic_parms = 0;
    ipopt_get_default_parameters(server,(SlvClientToken)sys,&(sys->p));
    sys->p.whose = (*statusindex);

    sys->presolved = 0;
    sys->resolve = 0;

    sys->n = -1;
    sys->m = -1;

    sys->s.ok = TRUE;
    sys->s.calc_ok = TRUE;
    sys->s.costsize = 0;
    sys->s.cost = NULL; /*redundant, but sanity preserving */

    sys->vlist = slv_get_solvers_var_list(server);
    sys->rlist = slv_get_solvers_rel_list(server);
    sys->obj = slv_get_obj_relation(server);
    if(sys->vlist == NULL) {
        ASC_FREE(sys);
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"IPOPT called with no variables.");
        *statusindex = -2;
        return NULL;
    }
    if(sys->rlist == NULL && sys->obj == NULL) {
        ASC_FREE(sys);
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"IPOPT called with no relations or objective.");
        *statusindex = -1;
        return NULL;
    }
    
    /* do nothing with the objective list, pars, bounds, extrels, etc etc */

    slv_check_var_initialization(server);
    *statusindex = 0;
    return((SlvClientToken)sys);
}

static int32 ipopt_destroy(slv_system_t server, SlvClientToken asys){
    UNUSED_PARAMETER(server);
    ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
    return 1;
}   

static int32 ipopt_eligible_solver(slv_system_t server){
    UNUSED_PARAMETER(server);
    ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
    return 1;
}

/*------------------------------------------------------------------------------
  SOLVER PARAMETERS
*/

static
int32 ipopt_get_default_parameters(slv_system_t server, SlvClientToken asys
        ,slv_parameters_t *parameters
){
    IpoptSystem *sys = NULL;
    struct slv_parameter *new_parms = NULL;
    int32 make_macros = 0;

    if(server != NULL && asys != NULL) {
        sys = SYS(asys);
        make_macros = 1;
    }

    if(parameters->parms == NULL) {
        new_parms = ASC_NEW_ARRAY_OR_NULL(struct slv_parameter,IPOPT_PARAMS);
        if(new_parms == NULL) {
            return -1;
        }
        parameters->parms = new_parms;
        parameters->dynamic_parms = 1;
    }

    parameters->num_parms = 0;

    slv_param_int(parameters,IPOPT_PARAM_MAX_ITER
        ,(SlvParameterInitInt){{"max_iter"
            ,"Maximum number of iterations",2
            ,"The algorithm terminates with an error message if the number of iterations exceeded this number."
        }, 3000, 0, 100000000}
    );

    slv_param_real(parameters,IPOPT_PARAM_TOL
        ,(SlvParameterInitReal){{"tol"
            ,"Desired convergence tolerance (relative)",2
            ,"Determines the convergence tolerance for the algorithm. The algorithm"
            " terminates successfully, if the (scaled) NLP error becomes smaller"
            " than this value, and if the (absolute) criteria according to "
            " 'dual_inf_tol', 'primal_inf_tol', and 'cmpl_inf_tol' are met. (This"
            " is epsilon_tol in Eqn. (6) in implementation paper). See also "
            " 'acceptable_tol' as a second termination criterion. Note, some other"
            " algorithmic features also use this quantity to determine thresholds"
            " etc."
        }, 1e-8, 0, 1e20}
    );

    slv_param_char(parameters,IPOPT_PARAM_MU_STRATEGY
        ,(SlvParameterInitChar){{"mu_strategy"
            ,"Update strategy for barrier parameter",6
            ,"Determines which barrier parameter update strategy is to be used."
            " 'MONOTONE' is the monotone (Fiacco-McCormick) strategy;"
            " 'ADAPTIVE' is the adaptive update strategy."
        }, "MONOTONE"}, (char *[]){
            "MONOTONE","ADAPTIVE",NULL
        }
    );

    asc_assert(parameters->num_parms==IPOPT_PARAMS);

    return 1;
}

static void ipopt_get_parameters(slv_system_t server, SlvClientToken asys
        , slv_parameters_t *parameters
){
    IpoptSystem *sys;
    UNUSED_PARAMETER(server);

    sys = SYS(asys);
    //if(check_system(sys)) return;
    mem_copy_cast(&(sys->p),parameters,sizeof(slv_parameters_t));
}


static void ipopt_set_parameters(slv_system_t server, SlvClientToken asys
        ,slv_parameters_t *parameters
){
    IpoptSystem *sys;
    UNUSED_PARAMETER(server);

    sys = SYS(asys);
    //if (check_system(sys)) return;
    mem_copy_cast(parameters,&(sys->p),sizeof(slv_parameters_t));
}

/*------------------------------------------------------------------------------
  EVALUATION FUNCTIONS
*/

/**
    update the model with new 'x' vector.
    @return 0 on success.
*/
int ipopt_update_model(IpoptSystem *sys, const double *x){
    ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
    return 1; /* error */
}

/** Function to evaluate the objective function f(x).
    @return 1 on success, 0 on failure

    @param n (in), the number of variables in the problem (dimension of 'x').
    @param x (in), the values for the primal variables, 'x' , at which 'f(x)' is to be evaluated.
    @param new_x (in), false if any evaluation method was previously called with the same values in 'x', true otherwise.
    @param obj_value (out) the value of the objective function ('f(x)').
*/
Bool ipopt_eval_f(Index n, Number *x, Bool new_x,  Number *obj_value, void *user_data){
    IpoptSystem *sys;
    sys = SYS(user_data);
    int res;

    asc_assert(n==sys->n);

    if(new_x){
        res = ipopt_update_model(sys,x);
        if(res)return 0; /* fail model update */
    }


    double val = 0;
    /* evaluate f(x) somehow */

    *obj_value = val;

    return 0; /* fail: not yet implemented */
}

Bool ipopt_eval_grad_f(Index n, Number* x, Bool new_x, Number* grad_f, void *user_data){
    IpoptSystem *sys;
    sys = SYS(user_data);
    int j, res;

    asc_assert(n==sys->n);

    if(new_x){
        res = ipopt_update_model(sys,x);
        if(res)return 0; /* fail model update */
    }


    /* evaluate grad_f(x) somehow */
    for(j=0; j<n; ++j){
        grad_f[j] = 0;
    }

    return 0; /* fail: not yet implemented */
}

Bool ipopt_eval_g(Index n, Number* x, Bool new_x, Index m, Number *g, void *user_data){
    IpoptSystem *sys;
    sys = SYS(user_data);
    int i, res;

    asc_assert(n==sys->n);
    asc_assert(m==sys->m);

    if(new_x){
        res = ipopt_update_model(sys,x);
        if(res)return 0; /* fail model update */
    }

    for(i=0; i<m; ++i){
        g[i] = 0;
    }

    return 0; /* fail: not yet implemented */
}

Bool ipopt_eval_jac_g(Index n, Number* x, Bool new_x, Index m
        , Index nele_jac, Index* iRow, Index *jCol, Number* values
        , void *user_data
){
    IpoptSystem *sys;
    sys = SYS(user_data);
    int i,j,res;

    asc_assert(n==sys->n);
    asc_assert(nele_jac==sys->nnzJ);
    asc_assert(m==sys->m);

    if(new_x){
        res = ipopt_update_model(sys,x);
        if(res)return 0; /* fail model update */
    }

    /* loop through the constraints */
    for(i=0; i<m; ++i){
        /* get derivatives for constraint i */
        /* insert the derivatives into the matrix in row i, columns j */
    }

    return 0; /* fail: not yet implemented */
}

Bool ipopt_eval_h(Index n, Number* x, Bool new_x
        , Number obj_factor, Index m, Number* lambda
        , Bool new_lambda, Index nele_hess, Index* iRow
        , Index* jCol, Number* values
){
    /* not sure about this one yet: do all the 2nd deriv things work for this? */
    return 0; /* fail: not yet implemented */
}

/*------------------------------------------------------------------------------
  SOLVE ROUTINES
*/

static int ipopt_solve(slv_system_t server, SlvClientToken asys){
    IpoptSystem *sys;
    UNUSED_PARAMETER(server);
    int i;

    sys = SYS(asys);

    /* set the number of variables and allocate space for the bounds */
    sys->x_L = ASC_NEW_ARRAY(Number,sys->n);
    sys->x_U = ASC_NEW_ARRAY(Number,sys->n);

    /* set the values for the variable bounds */

#if 0
// sample code for the C interface

int main(){

  /* set the number of constraints and allocate space for the bounds */
  m=2;
  g_L = (Number*)malloc(sizeof(Number)*m);
  g_U = (Number*)malloc(sizeof(Number)*m);
  /* 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);
  
  /* We can free the memory now - the values for the bounds have been
     copied internally in CreateIpoptProblem */
  free(x_L);
  free(x_U);
  free(g_L);
  free(g_U);

  /* set some options */
  AddIpoptNumOption(nlp, "tol", 1e-9);
  AddIpoptStrOption(nlp, "mu_strategy", "adaptive");

  /* allocate space for the initial point and set the values */
  x = (Number*)malloc(sizeof(Number)*n);
  x[0] = 1.0;
  x[1] = 5.0;
  x[2] = 5.0;
  x[3] = 1.0;

  /* allocate space to store the bound multipliers at the solution */
  mult_x_L = (Number*)malloc(sizeof(Number)*n);
  mult_x_U = (Number*)malloc(sizeof(Number)*n);

  /* solve the problem */
  status = IpoptSolve(nlp, x, NULL, &obj, NULL, mult_x_L, mult_x_U, NULL);

  if (status == Solve_Succeeded) {
    printf("\n\nSolution of the primal variables, x\n");
    for (i=0; i<n; i++) {
      printf("x[%d] = %e\n", i, x[i]); 
    }

    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);

  return 0;
}
#endif


    ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
    return 1;
}

static int ipopt_presolve(slv_system_t server, SlvClientToken asys){
    IpoptSystem *sys;
    struct var_variable **vp;
    struct rel_relation **rp;
    int cap, ind;
    int matrix_creation_needed = 1;
    int *cntvect, temp;

    CONSOLE_DEBUG("PRESOLVE");

    sys = SYS(asys);
    //iteration_begins(sys);
    //check_system(sys);

    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;
    }

#if 0
    if(sys->presolved > 0) { /* system has been presolved before */
        if(!conopt_dof_changed(sys) /*no changes in fixed or included flags*/
                && sys->p.partition == sys->J.old_partition
                && sys->obj == sys->old_obj
        ){
            matrix_creation_needed = 0;
            CONOPT_CONSOLE_DEBUG("YOU JUST AVOIDED MATRIX DESTRUCTION/CREATION");
        }
    }
#endif

#if 0
    // check all this...

    /* 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); */
        }

        sys->presolved = 1; /* full presolve recognized here */
        sys->resolve = 0;   /* initialize resolve flag */

        /* provide nnz for jacobian matrix of constraints */

        /* provide sparsity structure for jacobian */

        /* provide nnz for hessian matrix */

        /* provide sparsity structure for hessian matrix */

        sys->J.old_partition = sys->p.partition;
        sys->old_obj = sys->obj;

        slv_sort_rels_and_vars(server,&(sys->con.m),&(sys->con.n));
        CONOPT_CONSOLE_DEBUG("FOUND %d CONSTRAINTS AND %d VARS",sys->con.m,sys->con.n);
        if (sys->obj != NULL) {
            CONOPT_CONSOLE_DEBUG("ADDING OBJECT AS A ROW");
            sys->con.m++; /* treat objective as a row */
        }

        cntvect = ASC_NEW_ARRAY(int,COIDEF_Size());
        COIDEF_Ini(cntvect);
        sys->con.cntvect = cntvect;
        CONOPT_CONSOLE_DEBUG("NUMBER OF CONSTRAINTS = %d",sys->con.m);
        COIDEF_NumVar(cntvect, &(sys->con.n));
        COIDEF_NumCon(cntvect, &(sys->con.m));
        sys->con.nz = num_jacobian_nonzeros(sys, &(sys->con.maxrow));
        COIDEF_NumNZ(cntvect, &(sys->con.nz));
        COIDEF_NumNlNz(cntvect, &(sys->con.nz));

        sys->con.base = 0;
        COIDEF_Base(cntvect,&(sys->con.base));
        COIDEF_ErrLim(cntvect, &(DOMLIM));
        COIDEF_ItLim(cntvect, &(ITER_LIMIT));

        if(sys->obj!=NULL){
            sys->con.optdir = relman_obj_direction(sys->obj);
            sys->con.objcon = sys->con.m - 1; /* objective will be last row */
            CONOPT_CONSOLE_DEBUG("SETTING OBJECTIVE CONSTRAINT TO BE %d",sys->con.objcon);
        }else{
            sys->con.optdir = 0;
            sys->con.objcon = 0;
        }
        COIDEF_OptDir(cntvect, &(sys->con.optdir));
        COIDEF_ObjCon(cntvect, &(sys->con.objcon));

        temp = 0;
        COIDEF_StdOut(cntvect, &temp);

        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);

        int debugfv = 1;
        COIDEF_DebugFV(cntvect, &debugfv);

        destroy_vectors(sys);
        destroy_matrices(sys);
        create_matrices(server,sys);
        create_vectors(sys);

        sys->s.block.current_reordered_block = -2;
    }

    /* 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;

    if( matrix_creation_needed ) {
        destroy_array(sys->s.cost);
        sys->s.cost = create_zero_array(sys->s.costsize,struct slv_block_cost);
        for( ind = 0; ind < sys->s.costsize; ++ind ) {
            sys->s.cost[ind].reorder_method = -1;
        }
    } else {
        reset_cost(sys->s.cost,sys->s.costsize);
    }

    /* set to go to first unconverged block */
    sys->s.block.current_block = -1;
    sys->s.block.current_size = 0;
    sys->s.calc_ok = TRUE;
    sys->s.block.iteration = 0;
    sys->objective =  MAXDOUBLE/2000.0;

    update_status(sys);
    iteration_ends(sys);
    sys->s.cost[sys->s.block.number_of].time=sys->s.cpu_elapsed;
#endif

    return 0;
}

static int ipopt_iterate(slv_system_t server, SlvClientToken asys){
    UNUSED_PARAMETER(server);
    ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
    return 1;
}

static int ipopt_resolve(slv_system_t server, SlvClientToken asys){
    UNUSED_PARAMETER(server);

    /* if implementing this, use the 'warm start' thing in IPOPT */

    /* provide initial values of the 'multipliers' */

    ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
    return 1;
}

static const SlvFunctionsT ipopt_internals = {
    67
    ,"IPOPT"
    ,ipopt_create
    ,ipopt_destroy
    ,ipopt_eligible_solver
    ,ipopt_get_default_parameters
    ,ipopt_get_parameters
    ,ipopt_set_parameters
    ,NULL
    ,ipopt_solve
    ,ipopt_presolve
    ,ipopt_iterate
    ,ipopt_resolve
    ,NULL
    ,NULL
    ,NULL
};

int ipopt_register(void){
    ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to load IPOPT");
    return 1;
    /* return solver_register(&tron_internals); */
}
1	jpye	1528	/* 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	jpye	1530	The IPOPT solver is documented at http://projects.coin-or.org/Ipopt/
26	jpye	1528	//
27			ASCEND wrapper for IPOPT originally by John Pye, Jun 2007.
28			*/
29
30	jpye	1529	#include <utilities/config.h>
31
32			#ifndef WITH_IPOPT
33			# error "WITH_IPOPT must be defined in order to build this."
34			#endif
35
36	jpye	1528	#include <solver/solver.h>
37
38	jpye	1530	#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	jpye	1528	#include <ipopt/IpStdCInterface.h>
53
54			ASC_DLLSPEC SolverRegisterFn ipopt_register;
55
56	jpye	1530	/*------------------------------------------------------------------------------
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_SAFECALC
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			/*
86			Solver information
87			*/
88			int32 presolved; /* ? Has the system been presolved */
89			int32 resolve; /* ? Has the system been resolved */
90			slv_parameters_t p; /* Parameters */
91			slv_status_t s; /* Status (as of iteration end) */
92
93			int32 cap; /* Order of matrix/vectors */
94			int32 rank; /* Symbolic rank of problem */
95			int32 vused; /* Free and incident variables */
96			int32 vtot; /* length of varlist */
97			int32 rused; /* Included relations */
98			int32 rtot; /* length of rellist */
99			double clock; /* CPU time */
100
101			void *parm_array[IPOPT_PARAMS];
102			struct slv_parameter pa[IPOPT_PARAMS];
103
104			/*
105			IPOPT DATA
106			*/
107			Index n; /* number of variables */
108			Index m; /* number of constraints */
109
110			int nnzJ; /* number of non zeros in the jacobian of the constraints */
111
112			Number* x_L; /* lower bounds on x */
113			Number* x_U; /* upper bounds on x */
114			Number* g_L; /* lower bounds on g */
115			Number* g_U; /* upper bounds on g */
116
117			IpoptProblem nlp; /* IpoptProblem */
118
119			enum ApplicationReturnStatus status; /* Solve return code */
120			Number* x; /* starting point and solution vector */
121			Number* mult_x_L; /* lower bound multipliers at the solution */
122			Number* mult_x_U; /* upper bound multipliers at the solution */
123			Number objval; /* objective value */
124			Index i; /* generic counter */
125			};
126
127			typedef struct IpoptSystemStruct IpoptSystem;
128
129			static int ipopt_get_default_parameters(slv_system_t server, SlvClientToken asys, slv_parameters_t *parameters);
130
131			/*------------------------------------------------------------------------------
132			SYSTEM SETUP/DESTROY AND SOLVER ELIGIBILITY
133			*/
134
135			static SlvClientToken ipopt_create(slv_system_t server, int32*statusindex){
136			IpoptSystem *sys;
137
138			sys = ASC_NEW_CLEAR(IpoptSystem);
139			if(sys==NULL){
140			*statusindex = 1;
141			return sys;
142			}
143
144			sys->p.parms = sys->pa;
145			sys->p.dynamic_parms = 0;
146			ipopt_get_default_parameters(server,(SlvClientToken)sys,&(sys->p));
147			sys->p.whose = (*statusindex);
148
149			sys->presolved = 0;
150			sys->resolve = 0;
151
152			sys->n = -1;
153			sys->m = -1;
154
155			sys->s.ok = TRUE;
156			sys->s.calc_ok = TRUE;
157			sys->s.costsize = 0;
158			sys->s.cost = NULL; /redundant, but sanity preserving /
159
160			sys->vlist = slv_get_solvers_var_list(server);
161			sys->rlist = slv_get_solvers_rel_list(server);
162			sys->obj = slv_get_obj_relation(server);
163			if(sys->vlist == NULL) {
164			ASC_FREE(sys);
165			ERROR_REPORTER_HERE(ASC_PROG_ERR,"IPOPT called with no variables.");
166			*statusindex = -2;
167			return NULL;
168			}
169			if(sys->rlist == NULL && sys->obj == NULL) {
170			ASC_FREE(sys);
171			ERROR_REPORTER_HERE(ASC_PROG_ERR,"IPOPT called with no relations or objective.");
172			*statusindex = -1;
173			return NULL;
174			}
175
176			/* do nothing with the objective list, pars, bounds, extrels, etc etc */
177
178			slv_check_var_initialization(server);
179			*statusindex = 0;
180			return((SlvClientToken)sys);
181			}
182
183			static int32 ipopt_destroy(slv_system_t server, SlvClientToken asys){
184			UNUSED_PARAMETER(server);
185			ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
186			return 1;
187			}
188
189			static int32 ipopt_eligible_solver(slv_system_t server){
190			UNUSED_PARAMETER(server);
191			ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
192			return 1;
193			}
194
195			/*------------------------------------------------------------------------------
196			SOLVER PARAMETERS
197			*/
198
199			static
200			int32 ipopt_get_default_parameters(slv_system_t server, SlvClientToken asys
201			,slv_parameters_t *parameters
202			){
203			IpoptSystem *sys = NULL;
204			struct slv_parameter *new_parms = NULL;
205			int32 make_macros = 0;
206
207			if(server != NULL && asys != NULL) {
208			sys = SYS(asys);
209			make_macros = 1;
210			}
211
212			if(parameters->parms == NULL) {
213			new_parms = ASC_NEW_ARRAY_OR_NULL(struct slv_parameter,IPOPT_PARAMS);
214			if(new_parms == NULL) {
215			return -1;
216			}
217			parameters->parms = new_parms;
218			parameters->dynamic_parms = 1;
219			}
220
221			parameters->num_parms = 0;
222
223			slv_param_int(parameters,IPOPT_PARAM_MAX_ITER
224			,(SlvParameterInitInt){{"max_iter"
225			,"Maximum number of iterations",2
226			,"The algorithm terminates with an error message if the number of iterations exceeded this number."
227			}, 3000, 0, 100000000}
228			);
229
230			slv_param_real(parameters,IPOPT_PARAM_TOL
231			,(SlvParameterInitReal){{"tol"
232			,"Desired convergence tolerance (relative)",2
233			,"Determines the convergence tolerance for the algorithm. The algorithm"
234			" terminates successfully, if the (scaled) NLP error becomes smaller"
235			" than this value, and if the (absolute) criteria according to "
236			" 'dual_inf_tol', 'primal_inf_tol', and 'cmpl_inf_tol' are met. (This"
237			" is epsilon_tol in Eqn. (6) in implementation paper). See also "
238			" 'acceptable_tol' as a second termination criterion. Note, some other"
239			" algorithmic features also use this quantity to determine thresholds"
240			" etc."
241			}, 1e-8, 0, 1e20}
242			);
243
244			slv_param_char(parameters,IPOPT_PARAM_MU_STRATEGY
245			,(SlvParameterInitChar){{"mu_strategy"
246			,"Update strategy for barrier parameter",6
247			,"Determines which barrier parameter update strategy is to be used."
248			" 'MONOTONE' is the monotone (Fiacco-McCormick) strategy;"
249			" 'ADAPTIVE' is the adaptive update strategy."
250			}, "MONOTONE"}, (char *[]){
251			"MONOTONE","ADAPTIVE",NULL
252			}
253			);
254
255			asc_assert(parameters->num_parms==IPOPT_PARAMS);
256
257			return 1;
258			}
259
260			static void ipopt_get_parameters(slv_system_t server, SlvClientToken asys
261			, slv_parameters_t *parameters
262			){
263			IpoptSystem *sys;
264			UNUSED_PARAMETER(server);
265
266			sys = SYS(asys);
267			//if(check_system(sys)) return;
268			mem_copy_cast(&(sys->p),parameters,sizeof(slv_parameters_t));
269			}
270
271
272			static void ipopt_set_parameters(slv_system_t server, SlvClientToken asys
273			,slv_parameters_t *parameters
274			){
275			IpoptSystem *sys;
276			UNUSED_PARAMETER(server);
277
278			sys = SYS(asys);
279			//if (check_system(sys)) return;
280			mem_copy_cast(parameters,&(sys->p),sizeof(slv_parameters_t));
281			}
282
283			/*------------------------------------------------------------------------------
284			EVALUATION FUNCTIONS
285			*/
286
287			/**
288			update the model with new 'x' vector.
289			@return 0 on success.
290			*/
291			int ipopt_update_model(IpoptSystem sys, const double x){
292			ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
293			return 1; /* error */
294			}
295
296			/** Function to evaluate the objective function f(x).
297			@return 1 on success, 0 on failure
298
299			@param n (in), the number of variables in the problem (dimension of 'x').
300			@param x (in), the values for the primal variables, 'x' , at which 'f(x)' is to be evaluated.
301			@param new_x (in), false if any evaluation method was previously called with the same values in 'x', true otherwise.
302			@param obj_value (out) the value of the objective function ('f(x)').
303			*/
304			Bool ipopt_eval_f(Index n, Number x, Bool new_x, Number obj_value, void *user_data){
305			IpoptSystem *sys;
306			sys = SYS(user_data);
307			int res;
308
309			asc_assert(n==sys->n);
310
311			if(new_x){
312			res = ipopt_update_model(sys,x);
313			if(res)return 0; /* fail model update */
314			}
315
316
317			double val = 0;
318			/* evaluate f(x) somehow */
319
320			*obj_value = val;
321
322			return 0; /* fail: not yet implemented */
323			}
324
325			Bool ipopt_eval_grad_f(Index n, Number* x, Bool new_x, Number* grad_f, void *user_data){
326			IpoptSystem *sys;
327			sys = SYS(user_data);
328			int j, res;
329
330			asc_assert(n==sys->n);
331
332			if(new_x){
333			res = ipopt_update_model(sys,x);
334			if(res)return 0; /* fail model update */
335			}
336
337
338			/* evaluate grad_f(x) somehow */
339			for(j=0; j<n; ++j){
340			grad_f[j] = 0;
341			}
342
343			return 0; /* fail: not yet implemented */
344			}
345
346			Bool ipopt_eval_g(Index n, Number* x, Bool new_x, Index m, Number g, void user_data){
347			IpoptSystem *sys;
348			sys = SYS(user_data);
349			int i, res;
350
351			asc_assert(n==sys->n);
352			asc_assert(m==sys->m);
353
354			if(new_x){
355			res = ipopt_update_model(sys,x);
356			if(res)return 0; /* fail model update */
357			}
358
359			for(i=0; i<m; ++i){
360			g[i] = 0;
361			}
362
363			return 0; /* fail: not yet implemented */
364			}
365
366			Bool ipopt_eval_jac_g(Index n, Number* x, Bool new_x, Index m
367			, Index nele_jac, Index* iRow, Index jCol, Number values
368			, void *user_data
369			){
370			IpoptSystem *sys;
371			sys = SYS(user_data);
372			int i,j,res;
373
374			asc_assert(n==sys->n);
375			asc_assert(nele_jac==sys->nnzJ);
376			asc_assert(m==sys->m);
377
378			if(new_x){
379			res = ipopt_update_model(sys,x);
380			if(res)return 0; /* fail model update */
381			}
382
383			/* loop through the constraints */
384			for(i=0; i<m; ++i){
385			/* get derivatives for constraint i */
386			/* insert the derivatives into the matrix in row i, columns j */
387			}
388
389			return 0; /* fail: not yet implemented */
390			}
391
392			Bool ipopt_eval_h(Index n, Number* x, Bool new_x
393			, Number obj_factor, Index m, Number* lambda
394			, Bool new_lambda, Index nele_hess, Index* iRow
395			, Index* jCol, Number* values
396			){
397			/* not sure about this one yet: do all the 2nd deriv things work for this? */
398			return 0; /* fail: not yet implemented */
399			}
400
401			/*------------------------------------------------------------------------------
402			SOLVE ROUTINES
403			*/
404
405			static int ipopt_solve(slv_system_t server, SlvClientToken asys){
406			IpoptSystem *sys;
407			UNUSED_PARAMETER(server);
408			int i;
409
410			sys = SYS(asys);
411
412			/* set the number of variables and allocate space for the bounds */
413			sys->x_L = ASC_NEW_ARRAY(Number,sys->n);
414			sys->x_U = ASC_NEW_ARRAY(Number,sys->n);
415
416			/* set the values for the variable bounds */
417
418	jpye	1528	#if 0
419			// sample code for the C interface
420
421			int main(){
422
423			/* set the number of constraints and allocate space for the bounds */
424			m=2;
425			g_L = (Number)malloc(sizeof(Number)m);
426			g_U = (Number)malloc(sizeof(Number)m);
427			/* set the values of the constraint bounds */
428			g_L[0] = 25; g_U[0] = 2e19;
429			g_L[1] = 40; g_U[1] = 40;
430
431			/* create the IpoptProblem */
432			nlp = CreateIpoptProblem(n, x_L, x_U, m, g_L, g_U, 8, 10, 0,
433			&eval_f, &eval_g, &eval_grad_f,
434			&eval_jac_g, &eval_h);
435
436			/* We can free the memory now - the values for the bounds have been
437			copied internally in CreateIpoptProblem */
438			free(x_L);
439			free(x_U);
440			free(g_L);
441			free(g_U);
442
443			/* set some options */
444			AddIpoptNumOption(nlp, "tol", 1e-9);
445			AddIpoptStrOption(nlp, "mu_strategy", "adaptive");
446
447			/* allocate space for the initial point and set the values */
448			x = (Number)malloc(sizeof(Number)n);
449			x[0] = 1.0;
450			x[1] = 5.0;
451			x[2] = 5.0;
452			x[3] = 1.0;
453
454			/* allocate space to store the bound multipliers at the solution */
455			mult_x_L = (Number)malloc(sizeof(Number)n);
456			mult_x_U = (Number)malloc(sizeof(Number)n);
457
458			/* solve the problem */
459			status = IpoptSolve(nlp, x, NULL, &obj, NULL, mult_x_L, mult_x_U, NULL);
460
461			if (status == Solve_Succeeded) {
462			printf("\n\nSolution of the primal variables, x\n");
463			for (i=0; i<n; i++) {
464			printf("x[%d] = %e\n", i, x[i]);
465			}
466
467			printf("\n\nSolution of the bound multipliers, z_L and z_U\n");
468			for (i=0; i<n; i++) {
469			printf("z_L[%d] = %e\n", i, mult_x_L[i]);
470			}
471			for (i=0; i<n; i++) {
472			printf("z_U[%d] = %e\n", i, mult_x_U[i]);
473			}
474
475			printf("\n\nObjective value\n");
476			printf("f(x*) = %e\n", obj);
477			}
478
479			/* free allocated memory */
480			FreeIpoptProblem(nlp);
481			free(x);
482			free(mult_x_L);
483			free(mult_x_U);
484
485			return 0;
486			}
487			#endif
488
489	jpye	1530
490			ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
491			return 1;
492			}
493
494			static int ipopt_presolve(slv_system_t server, SlvClientToken asys){
495			IpoptSystem *sys;
496			struct var_variable **vp;
497			struct rel_relation **rp;
498			int cap, ind;
499			int matrix_creation_needed = 1;
500			int *cntvect, temp;
501
502			CONSOLE_DEBUG("PRESOLVE");
503
504			sys = SYS(asys);
505			//iteration_begins(sys);
506			//check_system(sys);
507
508			asc_assert(sys->vlist && sys->rlist);
509
510			sys->obj = slv_get_obj_relation(server); /may have changed objective/
511			if(!sys->obj){
512			ERROR_REPORTER_HERE(ASC_PROG_ERR,"No objective function was specified");
513			return -3;
514			}
515
516			#if 0
517			if(sys->presolved > 0) { /* system has been presolved before */
518			if(!conopt_dof_changed(sys) /no changes in fixed or included flags/
519			&& sys->p.partition == sys->J.old_partition
520			&& sys->obj == sys->old_obj
521			){
522			matrix_creation_needed = 0;
523			CONOPT_CONSOLE_DEBUG("YOU JUST AVOIDED MATRIX DESTRUCTION/CREATION");
524			}
525			}
526			#endif
527
528			#if 0
529			// check all this...
530
531			/* set all relations as being 'unsatisfied' to start with... */
532			rp=sys->rlist;
533			for( ind = 0; ind < sys->rtot; ++ind ) {
534			rel_set_satisfied(rp[ind],FALSE);
535			}
536
537			if( matrix_creation_needed ) {
538
539			cap = slv_get_num_solvers_rels(server);
540			sys->cap = slv_get_num_solvers_vars(server);
541			sys->cap = MAX(sys->cap,cap);
542			vp=sys->vlist;
543			for( ind = 0; ind < sys->vtot; ++ind ) {
544			var_set_in_block(vp[ind],FALSE);
545			}
546			rp=sys->rlist;
547			for( ind = 0; ind < sys->rtot; ++ind ) {
548			rel_set_in_block(rp[ind],FALSE);
549			/* rel_set_satisfied(rp[ind],FALSE); */
550			}
551
552			sys->presolved = 1; /* full presolve recognized here */
553			sys->resolve = 0; /* initialize resolve flag */
554
555			/* provide nnz for jacobian matrix of constraints */
556
557			/* provide sparsity structure for jacobian */
558
559			/* provide nnz for hessian matrix */
560
561			/* provide sparsity structure for hessian matrix */
562
563			sys->J.old_partition = sys->p.partition;
564			sys->old_obj = sys->obj;
565
566			slv_sort_rels_and_vars(server,&(sys->con.m),&(sys->con.n));
567			CONOPT_CONSOLE_DEBUG("FOUND %d CONSTRAINTS AND %d VARS",sys->con.m,sys->con.n);
568			if (sys->obj != NULL) {
569			CONOPT_CONSOLE_DEBUG("ADDING OBJECT AS A ROW");
570			sys->con.m++; /* treat objective as a row */
571			}
572
573			cntvect = ASC_NEW_ARRAY(int,COIDEF_Size());
574			COIDEF_Ini(cntvect);
575			sys->con.cntvect = cntvect;
576			CONOPT_CONSOLE_DEBUG("NUMBER OF CONSTRAINTS = %d",sys->con.m);
577			COIDEF_NumVar(cntvect, &(sys->con.n));
578			COIDEF_NumCon(cntvect, &(sys->con.m));
579			sys->con.nz = num_jacobian_nonzeros(sys, &(sys->con.maxrow));
580			COIDEF_NumNZ(cntvect, &(sys->con.nz));
581			COIDEF_NumNlNz(cntvect, &(sys->con.nz));
582
583			sys->con.base = 0;
584			COIDEF_Base(cntvect,&(sys->con.base));
585			COIDEF_ErrLim(cntvect, &(DOMLIM));
586			COIDEF_ItLim(cntvect, &(ITER_LIMIT));
587
588			if(sys->obj!=NULL){
589			sys->con.optdir = relman_obj_direction(sys->obj);
590			sys->con.objcon = sys->con.m - 1; /* objective will be last row */
591			CONOPT_CONSOLE_DEBUG("SETTING OBJECTIVE CONSTRAINT TO BE %d",sys->con.objcon);
592			}else{
593			sys->con.optdir = 0;
594			sys->con.objcon = 0;
595			}
596			COIDEF_OptDir(cntvect, &(sys->con.optdir));
597			COIDEF_ObjCon(cntvect, &(sys->con.objcon));
598
599			temp = 0;
600			COIDEF_StdOut(cntvect, &temp);
601
602			COIDEF_ReadMatrix(cntvect, &conopt_readmatrix);
603			COIDEF_FDEval(cntvect, &conopt_fdeval);
604			COIDEF_Option(cntvect, &conopt_option);
605			COIDEF_Solution(cntvect, &conopt_solution);
606			COIDEF_Status(cntvect, &conopt_status);
607			COIDEF_Message(cntvect, &asc_conopt_message);
608			COIDEF_ErrMsg(cntvect, &conopt_errmsg);
609			COIDEF_Progress(cntvect, &asc_conopt_progress);
610
611			int debugfv = 1;
612			COIDEF_DebugFV(cntvect, &debugfv);
613
614			destroy_vectors(sys);
615			destroy_matrices(sys);
616			create_matrices(server,sys);
617			create_vectors(sys);
618
619			sys->s.block.current_reordered_block = -2;
620			}
621
622			/* Reset status */
623			sys->con.optimized = 0;
624			sys->s.iteration = 0;
625			sys->s.cpu_elapsed = 0.0;
626			sys->s.converged = sys->s.diverged = sys->s.inconsistent = FALSE;
627			sys->s.block.previous_total_size = 0;
628			sys->s.costsize = 1+sys->s.block.number_of;
629
630			if( matrix_creation_needed ) {
631			destroy_array(sys->s.cost);
632			sys->s.cost = create_zero_array(sys->s.costsize,struct slv_block_cost);
633			for( ind = 0; ind < sys->s.costsize; ++ind ) {
634			sys->s.cost[ind].reorder_method = -1;
635			}
636			} else {
637			reset_cost(sys->s.cost,sys->s.costsize);
638			}
639
640			/* set to go to first unconverged block */
641			sys->s.block.current_block = -1;
642			sys->s.block.current_size = 0;
643			sys->s.calc_ok = TRUE;
644			sys->s.block.iteration = 0;
645			sys->objective = MAXDOUBLE/2000.0;
646
647			update_status(sys);
648			iteration_ends(sys);
649			sys->s.cost[sys->s.block.number_of].time=sys->s.cpu_elapsed;
650			#endif
651
652			return 0;
653			}
654
655			static int ipopt_iterate(slv_system_t server, SlvClientToken asys){
656			UNUSED_PARAMETER(server);
657			ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
658			return 1;
659			}
660
661			static int ipopt_resolve(slv_system_t server, SlvClientToken asys){
662			UNUSED_PARAMETER(server);
663
664			/* if implementing this, use the 'warm start' thing in IPOPT */
665
666			/* provide initial values of the 'multipliers' */
667
668			ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
669			return 1;
670			}
671
672			static const SlvFunctionsT ipopt_internals = {
673			67
674			,"IPOPT"
675			,ipopt_create
676			,ipopt_destroy
677			,ipopt_eligible_solver
678			,ipopt_get_default_parameters
679			,ipopt_get_parameters
680			,ipopt_set_parameters
681			,NULL
682			,ipopt_solve
683			,ipopt_presolve
684			,ipopt_iterate
685			,ipopt_resolve
686			,NULL
687			,NULL
688			,NULL
689			};
690
691	jpye	1528	int ipopt_register(void){
692			ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to load IPOPT");
693			return 1;
694			/* return solver_register(&tron_internals); */
695			}