From 4ebb38eee0e70bd709b56625124178d6be2d56bb Mon Sep 17 00:00:00 2001 From: Andrew Cron Date: Mon, 2 May 2016 21:48:37 -0400 Subject: [PATCH 1/2] added support for warmstart multipliers being passed to solver with example --- examples/hs071_PY3_warmstart.py | 166 ++++++++++++++++++++++++++++++++ setup.py | 3 +- src/pyipoptcoremodule.c | 138 ++++++++++++++++---------- 3 files changed, 257 insertions(+), 50 deletions(-) create mode 100644 examples/hs071_PY3_warmstart.py diff --git a/examples/hs071_PY3_warmstart.py b/examples/hs071_PY3_warmstart.py new file mode 100644 index 0000000..4b63f50 --- /dev/null +++ b/examples/hs071_PY3_warmstart.py @@ -0,0 +1,166 @@ +#!/usr/bin/python + +# Author: Eric Xu. Washington University +# The same model as Ipopt/examples/hs071 + +import pyipopt +from numpy import * + +def print_variable(variable_name, value): + for i in range(len(value)): + print("{} {}".format(variable_name + "["+str(i)+"] =", value[i])) + +nvar = 4 +x_L = ones((nvar), dtype=float_) * 1.0 +x_U = ones((nvar), dtype=float_) * 5.0 + +ncon = 2 + +g_L = array([25.0, 40.0]) +g_U = array([2.0*pow(10.0, 19), 40.0]) + +def eval_f(x, user_data = None): + assert len(x) == 4 + return x[0] * x[3] * (x[0] + x[1] + x[2]) + x[2] + +def eval_grad_f(x, user_data = None): + assert len(x) == 4 + grad_f = array([ + x[0] * x[3] + x[3] * (x[0] + x[1] + x[2]) , + x[0] * x[3], + x[0] * x[3] + 1.0, + x[0] * (x[0] + x[1] + x[2]) + ], float_) + return grad_f; + +def eval_g(x, user_data= None): + assert len(x) == 4 + return array([ + x[0] * x[1] * x[2] * x[3], + x[0]*x[0] + x[1]*x[1] + x[2]*x[2] + x[3]*x[3] + ], float_) + +nnzj = 8 +def eval_jac_g(x, flag, user_data = None): + if flag: + return (array([0, 0, 0, 0, 1, 1, 1, 1]), + array([0, 1, 2, 3, 0, 1, 2, 3])) + else: + assert len(x) == 4 + return array([ x[1]*x[2]*x[3], + x[0]*x[2]*x[3], + x[0]*x[1]*x[3], + x[0]*x[1]*x[2], + 2.0*x[0], + 2.0*x[1], + 2.0*x[2], + 2.0*x[3] ]) + +nnzh = 10 +def eval_h(x, lagrange, obj_factor, flag, user_data = None): + if flag: + hrow = [0, 1, 1, 2, 2, 2, 3, 3, 3, 3] + hcol = [0, 0, 1, 0, 1, 2, 0, 1, 2, 3] + return (array(hcol), array(hrow)) + else: + values = zeros((10), float_) + values[0] = obj_factor * (2*x[3]) + values[1] = obj_factor * (x[3]) + values[2] = 0 + values[3] = obj_factor * (x[3]) + values[4] = 0 + values[5] = 0 + values[6] = obj_factor * (2*x[0] + x[1] + x[2]) + values[7] = obj_factor * (x[0]) + values[8] = obj_factor * (x[0]) + values[9] = 0 + values[1] += lagrange[0] * (x[2] * x[3]) + + values[3] += lagrange[0] * (x[1] * x[3]) + values[4] += lagrange[0] * (x[0] * x[3]) + + values[6] += lagrange[0] * (x[1] * x[2]) + values[7] += lagrange[0] * (x[0] * x[2]) + values[8] += lagrange[0] * (x[0] * x[1]) + values[0] += lagrange[1] * 2 + values[2] += lagrange[1] * 2 + values[5] += lagrange[1] * 2 + values[9] += lagrange[1] * 2 + return values + +def apply_new(x): + return True + +nlp = pyipopt.create(nvar, x_L, x_U, ncon, g_L, g_U, nnzj, nnzh, eval_f, eval_grad_f, eval_g, eval_jac_g) + +x0 = array([1.0, 5.0, 5.0, 1.0]) +pi0 = array([1.0, 1.0]) + +""" +print x0 +print nvar, ncon, nnzj +print x_L, x_U +print g_L, g_U +print eval_f(x0) +print eval_grad_f(x0) +print eval_g(x0) +a = eval_jac_g(x0, True) +print "a = ", a[1], a[0] +print eval_jac_g(x0, False) +print eval_h(x0, pi0, 1.0, False) +print eval_h(x0, pi0, 1.0, True) +""" + +""" You can set Ipopt options by calling nlp.num_option, nlp.str_option +or nlp.int_option. For instance, to set the tolarance by calling + + nlp.num_option('tol', 1e-8) + +For a complete list of Ipopt options, refer to + + http://www.coin-or.org/Ipopt/documentation/node59.html + +Note that Ipopt distinguishs between Int, Num, and Str options, yet sometimes +does not explicitly tell you which option is which. If you are not sure about +the option's type, just try it in PyIpopt. If you try to set one type of +option using the wrong function, Pyipopt will remind you of it. """ + +print("Going to call solve for 4 iterations") +print("x0 = {}".format(x0)) +nlp.int_option('max_iter', 4) # limit the number of max iterations +x, zl, zu, constraint_multipliers, obj, status = nlp.solve(x0) +# import pdb; pdb.set_trace() +nlp.close() + +print("Solution of the bound multipliers, z_L and z_U") +print_variable("z_L", zl) +print_variable("z_U", zu) + +print("Solution of the constraint multipliers, lambda") +print_variable("lambda", constraint_multipliers) + + +nlp = pyipopt.create(nvar, x_L, x_U, ncon, g_L, g_U, nnzj, nnzh, eval_f, eval_grad_f, eval_g, eval_jac_g) +nlp.str_option('warm_start_init_point', 'yes') +nlp.num_option('warm_start_bound_push', 1e-8) +nlp.num_option('warm_start_slack_bound_push', 1e-8) +nlp.num_option('warm_start_mult_bound_push', 1e-8) +nlp.int_option('print_level', 5) +print("Starting at previous solution and solving again") +x, zl, zu, constraint_multipliers, obj, status = nlp.solve(x, mult_g=constraint_multipliers[:-1], + mult_x_L=zl, mult_x_U=zu) +nlp.close() + +print("Solution of the primal variables, x") +print_variable("x", x) + +print("Solution of the bound multipliers, z_L and z_U") +print_variable("z_L", zl) +print_variable("z_U", zu) + +print("Solution of the constraint multipliers, lambda") +print_variable("lambda", constraint_multipliers) + +print("Objective value") +print("f(x*) = {}".format(obj)) + diff --git a/setup.py b/setup.py index 5899b61..4c8c1e9 100644 --- a/setup.py +++ b/setup.py @@ -12,7 +12,8 @@ # for my linux distribution was buggy, # so by the time you read this the bugs have probably been fixed # and you will want to specify a different directory here. -IPOPT_DIR = '/usr/local/' +#IPOPT_DIR = '/usr/local/' +IPOPT_DIR = '/home/andrecr/store/local/ipopt' import os from distutils.core import setup diff --git a/src/pyipoptcoremodule.c b/src/pyipoptcoremodule.c index 8d37845..9cb6b86 100644 --- a/src/pyipoptcoremodule.c +++ b/src/pyipoptcoremodule.c @@ -16,12 +16,16 @@ * Let's put the static char docs at the beginning of this file... */ -static char PYIPOPT_SOLVE_DOC[] = "solve(x) -> (x, ml, mu, obj)\n \ +static char PYIPOPT_SOLVE_DOC[] = "solve(x, [mult_g, mult_x_L, mult_x_U]) -> (x, ml, mu, obj)\n \ \n \ Call Ipopt to solve problem created before and return \n \ a tuple that contains final solution x, upper and lower\n \ bound for multiplier, final objective function obj, \n \ - and the return status of ipopt. \n"; + and the return status of ipopt. \n \ + \n \ + mult_g, mult_x_L, mult_x_U are optional keyword only arguments \n \ + allowing previous values of bound multipliers to be passed in warm \n \ + start applications."; static char PYIPOPT_SET_INTERMEDIATE_CALLBACK_DOC[] = "set_intermediate_callback(callback_function)\n \ @@ -101,7 +105,7 @@ static void problem_dealloc(PyObject * self) Py_TYPE(self)->tp_free((PyObject*)self); } -PyObject *solve(PyObject * self, PyObject * args); +PyObject *solve(PyObject * self, PyObject * args, PyObject * keywds); PyObject *set_intermediate_callback(PyObject * self, PyObject * args); PyObject *close_model(PyObject * self, PyObject * args); @@ -174,7 +178,7 @@ static PyObject *add_num_option(PyObject * self, PyObject * args) } PyMethodDef problem_methods[] = { - {"solve", solve, METH_VARARGS, PYIPOPT_SOLVE_DOC} + {"solve", (PyCFunction)solve, METH_VARARGS | METH_KEYWORDS, PYIPOPT_SOLVE_DOC} , {"set_intermediate_callback", set_intermediate_callback, METH_VARARGS, PYIPOPT_SET_INTERMEDIATE_CALLBACK_DOC} @@ -551,8 +555,39 @@ PyObject *set_intermediate_callback(PyObject * self, PyObject * args) } } -PyObject *solve(PyObject * self, PyObject * args) +// too much cut and paste ... +#define SOLVE_CLEANUP() { \ + /* clean up and return */ \ + if (retval == NULL) { \ + Py_XDECREF(x); \ + Py_XDECREF(mL); \ + Py_XDECREF(mU); \ + Py_XDECREF(lambda); \ + } \ + SAFE_FREE(newx0); \ + return retval; \ + } + +#define SOLVE_CLEANUP_NULL() { \ + retval = NULL; \ + SOLVE_CLEANUP() \ +} + +#define SOLVE_CLEANUP_MEMORY() { \ + retval = PyErr_NoMemory(); \ + SOLVE_CLEANUP() \ +} + +#define SOLVE_CLEANUP_TYPE(err) { \ + PyErr_SetString(PyExc_TypeError, err); \ + retval = NULL; \ + SOLVE_CLEANUP() \ +} + +PyObject *solve(PyObject * self, PyObject * args, PyObject * keywds) { + static char *kwlist[] = {"x0", "userdata", "mult_g", "mult_x_L", "mult_x_U", NULL}; + enum ApplicationReturnStatus status; /* Solve return code */ int i; int n; @@ -569,6 +604,7 @@ PyObject *solve(PyObject * self, PyObject * args) npy_intp dlambda[1]; PyArrayObject *x = NULL, *mL = NULL, *mU = NULL, *lambda = NULL; + PyArrayObject *mL_in = NULL, *mU_in = NULL, *lambda_in = NULL; Number obj; /* objective value */ PyObject *retval = NULL; @@ -578,17 +614,12 @@ PyObject *solve(PyObject * self, PyObject * args) Number *newx0 = NULL; - if (!PyArg_ParseTuple(args, "O!|O", &PyArray_Type, &x0, &myuserdata)) { - retval = NULL; - /* clean up and return */ - if (retval == NULL) { - Py_XDECREF(x); - Py_XDECREF(mL); - Py_XDECREF(mU); - Py_XDECREF(lambda); - } - SAFE_FREE(newx0); - return retval; + if (!PyArg_ParseTupleAndKeywords(args, keywds, "O!|O$O!O!O!", kwlist, + &PyArray_Type, &x0, &myuserdata, + &PyArray_Type, &lambda_in, // mult_g + &PyArray_Type, &mL_in, // mult_x_L + &PyArray_Type, &mU_in)) { // mult_x_Y + SOLVE_CLEANUP_NULL() } if (myuserdata != NULL) { bigfield->userdata = myuserdata; @@ -598,18 +629,7 @@ PyObject *solve(PyObject * self, PyObject * args) */ } if (nlp == NULL) { - PyErr_SetString(PyExc_TypeError, - "nlp objective passed to solve is NULL\n Problem created?\n"); - retval = NULL; - /* clean up and return */ - if (retval == NULL) { - Py_XDECREF(x); - Py_XDECREF(mL); - Py_XDECREF(mU); - Py_XDECREF(lambda); - } - SAFE_FREE(newx0); - return retval; + SOLVE_CLEANUP_TYPE("nlp objective passed to solve is NULL\n Problem created?\n") } if (bigfield->eval_h_python == NULL) { AddIpoptStrOption(nlp, "hessian_approximation", "limited-memory"); @@ -622,41 +642,60 @@ PyObject *solve(PyObject * self, PyObject * args) x = (PyArrayObject *) PyArray_SimpleNew(1, dX, PyArray_DOUBLE); if (!x) { - retval = PyErr_NoMemory(); - /* clean up and return */ - if (retval == NULL) { - Py_XDECREF(x); - Py_XDECREF(mL); - Py_XDECREF(mU); - Py_XDECREF(lambda); - } - SAFE_FREE(newx0); - return retval; + SOLVE_CLEANUP_MEMORY() } newx0 = (Number *) malloc(sizeof(Number) * n); if (!newx0) { - retval = PyErr_NoMemory(); - /* clean up and return */ - if (retval == NULL) { - Py_XDECREF(x); - Py_XDECREF(mL); - Py_XDECREF(mU); - Py_XDECREF(lambda); - } - SAFE_FREE(newx0); - return retval; + SOLVE_CLEANUP_MEMORY() } double *xdata = (double *)x0->data; + double *ydata; for (i = 0; i < n; i++) newx0[i] = xdata[i]; /* Allocate multiplier arrays */ - + int num_passed = 0; // must pass all or none mL = (PyArrayObject *) PyArray_SimpleNew(1, dX, PyArray_DOUBLE); + if(mL_in != NULL){ + if(mL_in->dimensions[0] != n){ + SOLVE_CLEANUP_TYPE("mult_x_L must be the same length as x0.\n") + } + num_passed += 1; + xdata = (double *) mL->data; + ydata = (double *) mL_in->data; + for (i = 0; i < dX[0]; i++) + xdata[i] = ydata[i]; + } mU = (PyArrayObject *) PyArray_SimpleNew(1, dX, PyArray_DOUBLE); + if(mU_in != NULL){ + if(mU_in->dimensions[0] != n){ + SOLVE_CLEANUP_TYPE("mult_x_U must be the same length as x0.\n") + } + num_passed += 1; + xdata = (double *) mU->data; + ydata = (double *) mU_in->data; + for (i = 0; i < dX[0]; i++) + xdata[i] = ydata[i]; + } dlambda[0] = m; lambda = (PyArrayObject *) PyArray_SimpleNew(1, dlambda, PyArray_DOUBLE); + if(lambda_in != NULL){ + if(lambda_in->dimensions[0] != m){ + SOLVE_CLEANUP_TYPE("mult_g must be the same length as the constraints.\n") + } + num_passed += 1; + xdata = (double *) lambda->data; + ydata = (double *) lambda_in->data; + for (i = 0; i < dlambda[0]; i++) + xdata[i] = ydata[i]; + } + + // some error checking for warm start + if( (num_passed != 0) & (num_passed != 3) ){ + SOLVE_CLEANUP_TYPE("If passing multipliers, you must pass them all.\n") + } + /* For status code, see IpReturnCodes_inc.h in Ipopt */ @@ -682,6 +721,7 @@ PyObject *solve(PyObject * self, PyObject * args) Py_XDECREF(mU); Py_XDECREF(lambda); + SAFE_FREE(newx0); return retval; } From 51642eda63507988e8aba1bb7948160c5cec4677 Mon Sep 17 00:00:00 2001 From: Andrew Cron Date: Tue, 3 May 2016 13:31:04 -0400 Subject: [PATCH 2/2] Setup path --- setup.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/setup.py b/setup.py index 4c8c1e9..5899b61 100644 --- a/setup.py +++ b/setup.py @@ -12,8 +12,7 @@ # for my linux distribution was buggy, # so by the time you read this the bugs have probably been fixed # and you will want to specify a different directory here. -#IPOPT_DIR = '/usr/local/' -IPOPT_DIR = '/home/andrecr/store/local/ipopt' +IPOPT_DIR = '/usr/local/' import os from distutils.core import setup