Skip to content

Warm Start in NLP Solve #27

New issue

Have a question about this project? # for a free GitHub account to open an issue and contact its maintainers and the community.

#

By clicking “#”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? # to your account

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Warm Start in NLP Solve #27

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
4ebb38e
added support for warmstart multipliers being passed to solver with e…
May 3, 2016
51642ed
Setup path
andrewcron May 3, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
166 changes: 166 additions & 0 deletions examples/hs071_PY3_warmstart.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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')
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I apologize for getting back to this PR after, hmm, two years. 😞

Just merged another PR of python3 support, so we have one version of hs071 as example. I wonder if the following 5 lines can be added to example with an if toggle.

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

138 changes: 89 additions & 49 deletions src/pyipoptcoremodule.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 \
Expand Down Expand Up @@ -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);

Expand Down Expand Up @@ -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}
Expand Down Expand Up @@ -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() { \
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

thanks for making it a macro!

/* 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};
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see. there is a subtle API change here (from args to kwargs), let me actually play with it a bit before saying it is good or bad.


enum ApplicationReturnStatus status; /* Solve return code */
int i;
int n;
Expand All @@ -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;
Expand All @@ -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;
Expand All @@ -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");
Expand All @@ -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 */

Expand All @@ -682,6 +721,7 @@ PyObject *solve(PyObject * self, PyObject * args)
Py_XDECREF(mU);
Py_XDECREF(lambda);


SAFE_FREE(newx0);
return retval;
}
Expand Down