User-scaling with bounds relaxation

[bknueven]

Summary:

When using nlp_scaling_method=user-scaling, the variable bounds are relaxed by Eqn. (35) in the implementation paper and then scaled according to the provided factor. It would friendlier to some models I'm working with (and I think more natural) to do things in the reverse order: apply the variable scaling factor to the bounds, and then relax the bounds according to Eqn. (35).

I'm wondering if the current behavior is as-intended, in which case I will work around it, or if this is not as-intended and is a matter of a (hopefully) straightforward patch.

Details:

I'm attempting to use scaling factors to scale some variables in a large-scale NLP. I've reduced the issue down to a trivial linear program, expressed in Pyomo:

import pyomo.environ as pyo
m = pyo.ConcreteModel()
factor = 1.0e-12
m.x = pyo.Var(bounds=(0.5*factor, 1.5*factor), initialize=factor)
m.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT)
m.scaling_factor[m.x] = 1./factor
 
m.o = pyo.Objective(expr=m.x/factor)

ipopt = pyo.SolverFactory("ipopt")
ipopt.options["print_level"] = 8
ipopt.options["nlp_scaling_method"] = "user-scaling"
ipopt.solve(m, tee=True, symbolic_solver_labels=True)

print(f"optimal value: {m.x.value}")

Or in math:

minimize 1.0e+12*x  subject to x within [5.0e-13, 1.5e-12] 

Notice that I'm scaling the variable x by 1.0e+12 to attempt to bring its value close to 1.0, and utilizing the user-scaling nlp_scaling_method.

The optimal solution reported by Ipopt:

optimal value: -9.999499999997493e-09

which is well outside, in relative terms, to the original lower bound of 5.0e-13.

Looking at the log for Ipopt, what's apparently happening is that Eqn. (35) in the implementation paper is applied to the unscaled variable bounds, which are then scaled. Below is some sample output:

Ipopt 3.14.4: print_level=8
nlp_scaling_method=user-scaling


List of options:

                                    Name   Value                # times used
                      nlp_scaling_method = user-scaling              1
                             print_level = 8                         2

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

This is Ipopt version 3.14.4, running with linear solver MUMPS 5.2.1.

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        0

objective scaling factor = 1
x scaling provided
No c scaling provided
No d scaling provided
DenseVector "x scaling vector" with 1 elements:
x scaling vector[    1]{x}= 1.0000000000000000e+12
DenseVector "modified x_L scaled" with 1 elements:
modified x_L scaled[    1]{x}=-9.9995000000000000e+03
DenseVector "modified x_U scaled" with 1 elements:
modified x_U scaled[    1]{x}= 1.0001500000000000e+04
DenseVector "modified d_L scaled" with 0 elements:
DenseVector "modified d_U scaled" with 0 elements:
DenseVector "initial x unscaled" with 1 elements:
initial x unscaled[    1]{x}= 9.9999999999999998e-13
Initial values of x sufficiently inside the bounds.
Initial values of s sufficiently inside the bounds.
Total number of variables............................:        1
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        1
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

For this problem, it would be better to scale the variable bounds and then apply Eqn. (35) from the implementation paper. Doing this by-hand (setting bound_relax_factor=0 and applying Eqn. (35) manually to the model), I converge to an solution equivalent to being well-scaled in the first place: 4.999999925092976e-13.

[svigerske]

I cannot say whether things were intentional or not, but I presume that to change this, one would touch OrigIpoptNLP::InitializeStructures(). This is where the scaling is computed and the bounds are relaxed:

Ipopt/src/Algorithm/IpOrigIpoptNLP.cpp

Lines 231 to 348 in a62679d

	if( !warm_start_same_structure_ )
	{
	retValue = nlp_->GetSpaces(x_space_, c_space_, d_space_, x_l_space_, px_l_space_, x_u_space_, px_u_space_,
	d_l_space_, pd_l_space_, d_u_space_, pd_u_space_, jac_c_space_, jac_d_space_, h_space_);
	if( !retValue )
	{
	jnlst_->Printf(J_WARNING, J_INITIALIZATION, "GetSpaces method for the NLP returns false.\n");
	return false;
	}
	// Check if the Hessian space is actually a limited-memory
	// approximation. If so, get the required information from the
	// NLP and create an appropreate h_space
	if( hessian_approximation_ == LIMITED_MEMORY )
	{
	SmartPtr<VectorSpace> approx_vecspace;
	SmartPtr<Matrix> P_approx;
	if( hessian_approximation_space_ == NONLINEAR_VARS )
	{
	nlp_->GetQuasiNewtonApproximationSpaces(approx_vecspace, P_approx);
	}
	if( IsValid(approx_vecspace) )
	{
	DBG_ASSERT(IsValid(P_approx));
	h_space_ = new LowRankUpdateSymMatrixSpace(x_space_->Dim(), ConstPtr(P_approx), ConstPtr(approx_vecspace),
	true);
	jnlst_->Printf(J_DETAILED, J_INITIALIZATION,
	"Hessian approximation will be done in smaller space of dimension %" IPOPT_INDEX_FORMAT " (instead of %" IPOPT_INDEX_FORMAT ")\n\n",
	P_approx->NCols(), P_approx->NRows());
	}
	else
	{
	DBG_ASSERT(IsNull(P_approx));
	h_space_ = new LowRankUpdateSymMatrixSpace(x_space_->Dim(), ConstPtr(P_approx), ConstPtr(x_space_), true);
	jnlst_->Printf(J_DETAILED, J_INITIALIZATION,
	"Hessian approximation will be done in the space of all %" IPOPT_INDEX_FORMAT " x variables.\n\n", x_space_->Dim());
	}
	}
	// Create the bounds structures
	x_L = x_l_space_->MakeNew();
	Px_L = px_l_space_->MakeNew();
	x_U = x_u_space_->MakeNew();
	Px_U = px_u_space_->MakeNew();
	d_L = d_l_space_->MakeNew();
	Pd_L = pd_l_space_->MakeNew();
	d_U = d_u_space_->MakeNew();
	Pd_U = pd_u_space_->MakeNew();
	retValue = nlp_->GetBoundsInformation(*Px_L, *x_L, *Px_U, *x_U, *Pd_L, *d_L, *Pd_U, *d_U);
	if( !retValue )
	{
	return false;
	}
	NLP_scaling()->DetermineScaling(x_space_, c_space_, d_space_, jac_c_space_, jac_d_space_, h_space_,
	scaled_jac_c_space_, scaled_jac_d_space_, scaled_h_space_, *Px_L, *x_L, *Px_U, *x_U);
	if( x_space_->Dim() < c_space_->Dim() && x_space_->Dim() > 0 )
	{
	char msg[128];
	Snprintf(msg, 127, "Too few degrees of freedom: %" IPOPT_INDEX_FORMAT " equality constraints but only %" IPOPT_INDEX_FORMAT " variables",
	c_space_->Dim(), x_space_->Dim());
	THROW_EXCEPTION(TOO_FEW_DOF, msg);
	}
	// cannot have any null pointers, want zero length vectors
	// instead of null - this will later need to be changed for _h;
	DBG_DO(retValue = IsValid(x_space_) && IsValid(c_space_) && IsValid(d_space_) && IsValid(x_l_space_)
	&& IsValid(px_l_space_) && IsValid(x_u_space_) && IsValid(px_u_space_) && IsValid(d_u_space_)
	&& IsValid(pd_u_space_) && IsValid(d_l_space_) && IsValid(pd_l_space_) && IsValid(jac_c_space_)
	&& IsValid(jac_d_space_) && IsValid(h_space_) && IsValid(scaled_jac_c_space_) && IsValid(scaled_jac_d_space_)
	&& IsValid(scaled_h_space_));
	DBG_ASSERT(retValue && "Model cannot return null vector or matrix prototypes or spaces,"
	" please return zero length vectors instead");
	}
	else
	{
	ASSERT_EXCEPTION(IsValid(x_space_), INVALID_WARMSTART,
	"OrigIpoptNLP called with warm_start_same_structure, but the problem is solved for the first time.");
	// Create the bounds structures
	x_L = x_l_space_->MakeNew();
	Px_L = px_l_space_->MakeNew();
	x_U = x_u_space_->MakeNew();
	Px_U = px_u_space_->MakeNew();
	d_L = d_l_space_->MakeNew();
	Pd_L = pd_l_space_->MakeNew();
	d_U = d_u_space_->MakeNew();
	Pd_U = pd_u_space_->MakeNew();
	retValue = nlp_->GetBoundsInformation(*Px_L, *x_L, *Px_U, *x_U, *Pd_L, *d_L, *Pd_U, *d_U);
	if( !retValue )
	{
	return false;
	}
	}
	x_L->Print(*jnlst_, J_MOREVECTOR, J_INITIALIZATION, "original x_L unscaled");
	x_U->Print(*jnlst_, J_MOREVECTOR, J_INITIALIZATION, "original x_U unscaled");
	d_L->Print(*jnlst_, J_MOREVECTOR, J_INITIALIZATION, "original d_L unscaled");
	d_U->Print(*jnlst_, J_MOREVECTOR, J_INITIALIZATION, "original d_U unscaled");
	if( bound_relax_factor_ > 0.0 )
	{
	SmartPtr<Vector> tmp;
	tmp = x_L->MakeNewCopy();
	orig_x_L_ = ConstPtr(tmp);
	tmp = x_U->MakeNewCopy();
	orig_x_U_ = ConstPtr(tmp);
	relax_bounds(-bound_relax_factor_, *x_L);
	relax_bounds(bound_relax_factor_, *x_U);
	relax_bounds(-bound_relax_factor_, *d_L);
	relax_bounds(bound_relax_factor_, *d_U);
	}

To avoid too many changes, a possible way would be to compute the amount of bound relaxation on the scaled bounds and then scale them back to the original problem and apply them to the original bounds. There should be an option to choose between computing the bound relaxation on the original and the scaled bounds.