[Ipopt] Very slow convergence of IPopt
Jan gekheid
jangekheid at gmail.com
Tue Dec 1 10:35:46 EST 2009
Hello Andreas,
I'm using the Java interface of IPopt and I've implemented an automatic
differentiation library (JEP Java Expression Parser) to automatically do the
differentiation of the objective function and the constraints. I implemented
this in the given HS071 example so that I can check the results. All the
automatically differentiated formulas come out the same as in the original
problem. The model and the results are attached to this mail.
It gives the right results, but the convergence is very slow. It takes more
than 3000 iterations and even then it stops because the maximum number of
iterations is exceeded. It is exactly the same problem as the original (as
for the formulas and starting conditions), except that the values for the
jacobian etc are calculated with the automatically differentiated functions.
This implementation is just to see whether the combination of using IPopt en
JEP is possible (which it is apparently). But I want to implement a much
larger problem with very much larger formulas for the hessian entries. It
would take days to compute just one problem if it takes this much time. Is
there anything I can do about that?
With kind regards,
Jan
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://list.coin-or.org/pipermail/ipopt/attachments/20091201/aaa17c85/attachment-0001.html
-------------- next part --------------
package jIpOpt;
import org.coinor.Ipopt;
import org.lsmp.djep.djep.DJep;
import org.nfunk.jep.*;
public class JIpOptMain
{
public static void main(String[]args)
{
new Data();
new HS071();
}
}
class HS071 extends Ipopt
{
// Problem sizes
int n, m, nele_jac, nele_hess;
double[] x;
boolean DJEP = true;
boolean print = false;
DJep djep;
private double objFunVal;
private double[] gVal;
private double[] gradFVal;
private double[][] jacGVal;
private double[][] hesVal;
private double[][][] hesGVal;
private Node objFunNode;
private Node[] xNodes;
private Node[] gradFNodes;
private Node[] gNodes;
private Node[][] hesNodes;
private Node[][] jacGNodes;
private Node[][][] hesGNodes;
/**
* Initialize the bounds and create the native Ipopt problem.
*/
public HS071()
{
x = getInitialGuess();
djep = new DJep();
djep.addStandardFunctions();
djep.addStandardConstants();
djep.addComplex();
djep.setAllowUndeclared(true);
djep.setAllowAssignment(true);
djep.setImplicitMul(true);
djep.addStandardDiffRules();
/* Number of nonzeros in the Jacobian of the constraints */
nele_jac = 8;
/* Number of nonzeros in the Hessian of the Lagrangian (lower or
* upper triangual part only) */
nele_hess = 10;
/* set the number of variables and allocate space for the bounds */
n=4;
double x_L[] = {1,1,1,1};
double x_U[] = {5,5,5,5};
/* set the number of constraints and allocate space for the bounds */
m=2;
/* set the values of the constraint bounds */
double g_L[] = {25,40};
double g_U[] = {2e19,40};
/* Index style for the irow/jcol elements */
int index_style = Ipopt.C_STYLE;
gVal = new double[m];
gradFVal = new double[n];
jacGVal = new double[m][n];
hesVal = new double[n][n];
hesGVal = new double[m][n][n];
xNodes = new Node[n];
gradFNodes = new Node[n];
gNodes = new Node[m];
jacGNodes = new Node[m][n];
hesNodes = new Node[n][n];
hesGNodes = new Node[m][n][n];
if(DJEP)
createNodes();
/* create the IpoptProblem */
create(n, x_L, x_U, m, g_L, g_U, nele_jac, nele_hess, index_style);
solve(x);
}
public double[] getInitialGuess(){
/* allocate space for the initial point and set the values */
double x[] = new double[4];
x[0] = 1.0;
x[1] = 5.0;
x[2] = 5.0;
x[3] = 1.0;
return x;
}
protected boolean eval_f(int n, double[] x, boolean new_x, double[] obj_value) {
if(DJEP){
evalNodes(x);
obj_value[0] = objFunVal;
}
else
obj_value[0] = x[0] * x[3] * (x[0] + x[1] + x[2]) + x[2];
Data.output.println("Obj fun val = " + obj_value[0]);
return true;
}
protected boolean eval_grad_f(int n, double[] x, boolean new_x, double[] grad_f) {
Data.output.println("grad_f");
if(DJEP)
{
evalNodes(x);
for(int i = 0; i < n; i++)
grad_f[i] = gradFVal[i];
}
else
{
grad_f[0] = x[0] * x[3] + x[3] * (x[0] + x[1] + x[2]);
grad_f[1] = x[0] * x[3];
grad_f[2] = x[0] * x[3] + 1;
grad_f[3] = x[0] * (x[0] + x[1] + x[2]);
}
for(int i = 0; i < 4; i++)
Data.output.println(grad_f[i]);
Data.output.println();
return true;
}
protected boolean eval_g(int n, double[] x, boolean new_x, int m, double[] g) {
Data.output.println("g");
if(DJEP == true)
{
evalNodes(x);
for(int i = 0; i < m; i ++)
g[i] = gVal[i];
}
else
{
g[0] = x[0] * x[1] * x[2] * x[3];
g[1] = x[0]*x[0] + x[1]*x[1] + x[2]*x[2] + x[3]*x[3];
}
Data.output.println(g[0]);
Data.output.println(g[1]);
Data.output.println();
return true;
}
protected boolean eval_jac_g(int n, double[] x, boolean new_x,
int m, int nele_jac, int[] iRow, int[] jCol, double[] values) {
Data.output.println("jac_g");
int iter = 0;
if (values == null) {
/* return the structure of the jacobian */
/* this particular jacobian is dense */
iRow[0] = 0; jCol[0] = 0;
iRow[1] = 0; jCol[1] = 1;
iRow[2] = 0; jCol[2] = 2;
iRow[3] = 0; jCol[3] = 3;
iRow[4] = 1; jCol[4] = 0;
iRow[5] = 1; jCol[5] = 1;
iRow[6] = 1; jCol[6] = 2;
iRow[7] = 1; jCol[7] = 3;
}
else {
/* return the values of the jacobian of the constraints */
if(DJEP == true)
{
evalNodes(x);
for(int j = 0; j < m; j++)
for(int i = 0; i < n; i++)
values[iter++] = jacGVal[j][i];
}
else
{
values[0] = x[1]*x[2]*x[3]; /* 0,0 */
values[1] = x[0]*x[2]*x[3]; /* 0,1 */
values[2] = x[0]*x[1]*x[3]; /* 0,2 */
values[3] = x[0]*x[1]*x[2]; /* 0,3 */
values[4] = 2*x[0]; /* 1,0 */
values[5] = 2*x[1]; /* 1,1 */
values[6] = 2*x[2]; /* 1,2 */
values[7] = 2*x[3]; /* 1,3 */
}
for(int i = 0; i < 8; i++)
Data.output.println(values[i]);
}
Data.output.println();
return true;
}
protected boolean eval_h(int n, double[] x, boolean new_x, double obj_factor, int m, double[] lambda, boolean new_lambda, int nele_hess, int[] iRow, int[] jCol, double[] values) {
int idx = 0; /* nonzero element counter */
int row = 0; /* row counter for loop */
int col = 0; /* col counter for loop */
int iter = 0;
Data.output.println("hessian");
if (values == null) {
/* return the structure. This is a symmetric matrix, fill the lower left
* triangle only. */
/* the hessian for this problem is actually dense */
idx=0;
for (row = 0; row < n; row++) {
for (col = 0; col <= row; col++) {
iRow[idx] = row;
jCol[idx] = col;
idx++;
}
}
}
else {
if(DJEP == true)
{
evalNodes(x);
for(int i = 0; i < n; i++){
for(int j = 0; j <= i; j++){
values[iter] += obj_factor*hesVal[i][j];
for(int k = 0; k < m; k++)
values[iter] += lambda[k]*hesGVal[k][i][j];
iter++;
}
}
}
else
{
/* return the values. This is a symmetric matrix, fill the lower left
* triangle only */
/* fill the objective portion */
values[0] = obj_factor * (2*x[3]); /* 0,0 */
values[1] = obj_factor * (x[3]); /* 1,0 */
values[2] = 0; /* 1,1 */
values[3] = obj_factor * (x[3]); /* 2,0 */
values[4] = 0; /* 2,1 */
values[5] = 0; /* 2,2 */
values[6] = obj_factor * (2*x[0] + x[1] + x[2]); /* 3,0 */
values[7] = obj_factor * (x[0]); /* 3,1 */
values[8] = obj_factor * (x[0]); /* 3,2 */
values[9] = 0; /* 3,3 */
/* add the portion for the first constraint */
values[1] += lambda[0] * (x[2] * x[3]); /* 1,0 */
values[3] += lambda[0] * (x[1] * x[3]); /* 2,0 */
values[4] += lambda[0] * (x[0] * x[3]); /* 2,1 */
values[6] += lambda[0] * (x[1] * x[2]); /* 3,0 */
values[7] += lambda[0] * (x[0] * x[2]); /* 3,1 */
values[8] += lambda[0] * (x[0] * x[1]); /* 3,2 */
/* add the portion for the second constraint */
values[0] += lambda[1] * 2; /* 0,0 */
values[2] += lambda[1] * 2; /* 1,1 */
values[5] += lambda[1] * 2; /* 2,2 */
values[9] += lambda[1] * 2; /* 3,3 */
for(int i = 0; i < 10; i++)
Data.output.println(values[i]);
}
}
Data.output.println();
return true;
}
public void createNodes()
{
double x[] = getInitialGuess();
try {
//determine xNodes
for(int i = 0; i < n; i++){
xNodes[i] = djep.preprocess(djep.parse("x"+i+"="+x[i]));
if(print)
djep.println(xNodes[i], Data.output);
}
//determine value function node
objFunNode = djep.preprocess(djep.parse("F=x0*x3*(x0+x1+x2)+x2"));
if(print){
djep.println(objFunNode, Data.output);
Data.output.println();
}
//determine gradF nodes
for(int i = 0; i < n; i++){
gradFNodes[i] = djep.simplify(djep.preprocess(djep.parse("diff(F,x"+i+")")));
if(print)
djep.println(gradFNodes[i], Data.output);
}
//determine gNodes
gNodes[0] = djep.preprocess(djep.parse("g0=x0*x1*x2*x3"));
gNodes[1] = djep.preprocess(djep.parse("g1=x0*x0+x1*x1+x2*x2+x3*x3"));
if(print){
Data.output.println();
djep.println(gNodes[0], Data.output);
djep.println(gNodes[1], Data.output);
Data.output.println();
}
//determine jacGNodes
for(int j = 0; j < m; j++){
for(int i = 0; i < n; i++){
jacGNodes[j][i] = djep.simplify(djep.preprocess(djep.parse("diff(g"+j+",x"+i+")")));
if(print){
djep.println(jacGNodes[j][i], Data.output);
}
}
}
if(print)
Data.output.println();
if(print)
Data.output.println("hesGNodes");
//determine hesGNodes
for(int k = 0; k < m; k++){
for(int i = 0; i < n; i++){
for(int j = 0;j < n; j++){
hesGNodes[k][i][j] = djep.simplify(djep.preprocess(djep.differentiate(jacGNodes[k][i], "x"+j)));
if(print){
djep.print(hesGNodes[k][i][j],Data.output);
Data.output.print("\t");
}
}
if(print)
Data.output.println();
}
if(print)
Data.output.println();
}
if(print)
Data.output.println("hesNodes");
//determine hesNodes
for(int i = 0; i < n; i++){
for(int j = 0; j < n; j++){
hesNodes[i][j] = djep.simplify(djep.preprocess(djep.parse("diff(diff(F"+",x"+i+"),x"+j+")")));
if(print){
djep.print(hesNodes[i][j],Data.output);
Data.output.print("\t\t\t");
}
}
if(print)
Data.output.println();
}
}
catch (ParseException e) {
e.printStackTrace();
}
}
public void evalNodes(double x[])
{
for(int i = 0; i < n; i++)
djep.setVarValue("x"+i, x[i]);
try{
objFunVal = (Double)djep.evaluate(objFunNode);
for(int i = 0; i < n; i++)
gradFVal[i] = (Double)djep.evaluate(gradFNodes[i]);
for(int i = 0; i < m; i++)
gVal[i] = (Double)djep.evaluate(gNodes[i]);
for(int j = 0; j < m; j++)
for(int i = 0; i < n; i++)
jacGVal[j][i] = (Double)djep.evaluate(jacGNodes[j][i]);
for(int k = 0; k < m; k++)
for(int i = 0; i < n; i++)
for(int j = 0;j < n; j++)
hesGVal[k][i][j] = (Double)djep.evaluate(hesGNodes[k][i][j]);
//determine hesNodes
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
hesVal[i][j] = (Double)djep.evaluate(hesNodes[i][j]);
}
catch (ParseException e) {
e.printStackTrace();
}
}
}
/* Original problem
* ******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Common Public License (CPL).
For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************
Number of nonzeros in equality constraint Jacobian...: 4
Number of nonzeros in inequality constraint Jacobian.: 4
Number of nonzeros in Lagrangian Hessian.............: 10
Total number of variables............................: 4
variables with only lower bounds: 0
variables with lower and upper bounds: 4
variables with only upper bounds: 0
Total number of equality constraints.................: 1
Total number of inequality constraints...............: 1
inequality constraints with only lower bounds: 1
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.6109693e+001 1.12e+001 5.28e-001 -1.0 0.00e+000 - 0.00e+000 0.00e+000 0
1 1.6982239e+001 7.30e-001 1.02e+001 -1.0 6.11e-001 - 7.19e-002 1.00e+000f 1
2 1.7318411e+001 6.94e-002 5.05e-001 -1.0 1.61e-001 - 1.00e+000 1.00e+000h 1
3 1.6849424e+001 3.15e-001 6.68e-002 -1.7 2.85e-001 - 7.94e-001 1.00e+000h 1
4 1.7051199e+001 4.79e-003 2.78e-003 -1.7 6.06e-002 - 1.00e+000 1.00e+000h 1
5 1.7011979e+001 7.20e-003 8.50e-003 -3.8 3.66e-002 - 9.45e-001 9.98e-001h 1
6 1.7014271e+001 1.16e-004 9.78e-006 -3.8 3.33e-003 - 1.00e+000 1.00e+000h 1
7 1.7014021e+001 7.04e-007 1.82e-007 -5.7 2.69e-004 - 1.00e+000 1.00e+000h 1
8 1.7014017e+001 9.82e-011 2.52e-011 -8.6 3.32e-006 - 1.00e+000 1.00e+000h 1
Number of Iterations....: 8
(scaled) (unscaled)
Objective...............: 1.7014017145179164e+001 1.7014017145179164e+001
Dual infeasibility......: 2.5164712419785811e-011 2.5164712419785811e-011
Constraint violation....: 1.7713830402499298e-011 1.7713830402499298e-011
Complementarity.........: 2.5277100427932871e-009 2.5277100427932871e-009
Overall NLP error.......: 2.5277100427932871e-009 2.5277100427932871e-009
Number of objective function evaluations = 9
Number of objective gradient evaluations = 9
Number of equality constraint evaluations = 9
Number of inequality constraint evaluations = 9
Number of equality constraint Jacobian evaluations = 9
Number of inequality constraint Jacobian evaluations = 9
Number of Lagrangian Hessian evaluations = 8
Total CPU secs in IPOPT (w/o function evaluations) = 0.204
Total CPU secs in NLP function evaluations = 0.015
EXIT: Optimal Solution Found.
%%%%%%%%%%%%%%%%%%%%% JEP implemented %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
3000 1.7014021e+001 7.82e-014 7.17e-004 -5.7 2.11e-007 - 1.00e+000 1.00e+000h 1
Number of Iterations....: 3000
(scaled) (unscaled)
Objective...............: 1.7014021213393679e+001 1.7014021213393679e+001
Dual infeasibility......: 7.1739878140864022e-004 7.1739878140864022e-004
Constraint violation....: 7.8159700933611020e-014 7.8159700933611020e-014
Complementarity.........: 1.8449144633328430e-006 1.8449144633328430e-006
Overall NLP error.......: 7.1739878140864022e-004 7.1739878140864022e-004
Number of objective function evaluations = 3001
Number of objective gradient evaluations = 3001
Number of equality constraint evaluations = 3001
Number of inequality constraint evaluations = 3001
Number of equality constraint Jacobian evaluations = 3001
Number of inequality constraint Jacobian evaluations = 3001
Number of Lagrangian Hessian evaluations = 3000
Total CPU secs in IPOPT (w/o function evaluations) = 6.389
Total CPU secs in NLP function evaluations = 7.517
EXIT: Maximum Number of Iterations Exceeded.
*/
More information about the Ipopt
mailing list