Best Practices for checking that AMPL found a valid solution

[brian]

I'm calling AMPL from my program and I wonder what is the best
practice to check that AMPL found a valid solution.

To get the optimization variables I use the <display> command to write
the vars to the log file. After AMPL is done I parse the log file.

So far to check that AMPL found a solution I've came up with:
1. Check AMPL process exit code is 0
2. In log file check that <objective>.result = solved

Are these checks good enough to test if AMPL found a valid solution?

I'm using windows vista and plan to eventually run on Red Hat.

Brian

[Paul]

You didn't specify LP v. IP/MIP v. NLP. By and large, I think what
you've got should suffice. If you are worried about possible
tolerance/premature convergence issues, you might display slacks that
have the wrong sign (or are nonzero in equation constraints), or
display the final MIP gap (hopefully zero, but most solvers have a
nonzero default optimality tolerance). There's another thread
somewhere around here about getting the final gap from CPLEX. I'm not
sure if CPLEX (or Gurobi) gives you the gap directly -- you may need
to get the final bound and subtract it from the "optimal" objective
value.

/Paul

[brian]

Hi Paul,

The problem is a NLP. I'm currently using minos for the solver and
plan on using multiple NLP solvers. I didn't specify the type of
problem because I assumed the same checks would work for all the types
of problems.

To get the results back into my program at the end of the AMPL model I
would do a bunch of displays for each var. In the process of creating
the model there were times I would make a mistake and the optimization
would not run but the display would dump out the original vals. I'm
trying to be careful so that I don't get an invalid solution and
interpret it as a valid solution.

I originally started the research as a masters student in electrical
engineering and I'm continuing the research as a PhD student. Below
is the model I used for my masters. I've wrapped the display vars
parts with landmarks to make parsing the output file easier. The
model will change a lot for my PhD.

Thanks for your help!!!!

Brian

----------------------------------------------------------------------------
option solver minos;

#no constants

param LdNmos := 4e-08;
param LdPmos := 5e-08;
param minWidth := 4e-07;
param minLength := 3.5e-07;
param maxWidth := 0.0002;
param maxLength := 2e-05;
param out := 1.65;
param MpCurFac := 1.144981598e-05;
param MpVth := 0.6798;
param MnCurFac := 6.455286683e-05;
param MnVth := 0.5591;
param Vnbias := 1;
param Vdd := 3.3;

var MpW:= 1e-06;
var MpL:= 1e-06;
var MnW:= 1e-06;
var MnL:= 1e-06;
var Vpbias:= 1.353;
var MnId= MnCurFac*(MnW/(MnL-2*LdNmos))*(Vnbias-MnVth)^2;
var MpId= MpCurFac*(MpW/(MpL-2*LdPmos))*(Vdd-Vpbias-MpVth)^2;

s.t. MinMpW : MpW >= minWidth;
s.t. MaxMpW : MpW <= maxWidth;
s.t. MinMpL : MpL >= minLength;
s.t. MaxMpL : MpL <= maxLength;
s.t. MinMnW : MnW >= minWidth;
s.t. MaxMnW : MnW <= maxWidth;
s.t. MinMnL : MnL >= minLength;
s.t. MaxMnL : MnL <= maxLength;
s.t. SatP : out <= Vpbias + MpVth;
s.t. SatP2 : Vdd-Vpbias >= MpVth;
s.t. EqualCurrent : MnId == MpId;

minimize area: MnW*MnL+MpW*MpL;

solve;

printf "startWriteVars\n";
display MpW;
display MpL;
display MnW;
display MnL;
display Vpbias;
display MnId;
display MpId;
printf "endWriteVars\n";

display area.result;
display area.result_num;
display area.message;
display area.exitcode;
------------------------------------------------------

[Paul]

Brian,

On Nov 14, 9:48 pm, brian <brian.down...@yakkadesign.com> wrote:
> Hi Paul,
>
> The problem is a NLP.  I'm currently using minos for the solver and
> plan on using multiple NLP solvers.  I didn't specify the type of
> problem because I assumed the same checks would work for all the types
> of problems.

Some apply across the board, some don't. For instance, with a MILP
you get a gap value that tells you how far your objective was from the
best bound (assuming you don't exhaust the search tree within your run
limits). If you do exhaust the search tree, then (barring numerical
instability problems) you're pretty much guaranteed to have an
optimum. With an NLP, in contrast, you always hit some sort of limit
(time, objective tolerance, whatever), and you need to pay more
attention to whether you got a genuine optimum. (And if your problem
is nonconvex -- which yours is -- then you really cannot be sure you
converged to a global optimum. See below.)

>
> To get the results back into my program at the end of the AMPL model I
> would do a bunch of displays for each var.  In the process of creating
> the model there were times I would make a mistake and the optimization
> would not run but the display would dump out the original vals.  I'm
> trying to be careful so that I don't get an invalid solution and
> interpret it as a valid solution.

To avoid getting the starting values, you could just check (after the
solve) that the terminal value for at least one variable is different
(by more than rounding errors) from the starting value.

For grins, I ran your model twice, the first time with the starting
values you specified. That run produced an "optimal" area of
2.284730862e-13. Then I shifted the starting values of MpW and MnW to
maxWidth (0.0002) and the starting values of MpL and MnL to maxLength
(2e-05) and reran the model. That run produced an "optimal" area of
1.031603201e-11 (almost 50 times larger!). Both runs produced the
"solved" result.

I don't think that MnId and MpId are convex functions of the other
variables, but I'm too lazy to check, particularly since the
EqualCurrent constraint (which sets two nonlinear functions equal to
each other) pretty much guarantees the feasible region will be
nonconvex. (It will be convex only if either lots of stuff magically
cancels leaving the remaining expressions linear or the feasible
region reduces to a single point.)

There are solvers out there for nonconvex nonlinear problems, but I
don't know which (if any) AMPL supports, since I avoid anything that
weird. I think global optimization requires multiple restarts of the
problem, and AFAIK there is nothing remotely resembling a guarantee of
optimality (or even near optimality). The best I can say is that,
unless I'm really off track here, you can't count on MINOS to give you
a solution you can take to the bank.

/Paul

[Michael Saunders]

Brian and Paul,

Someone should mention the importance of scaling here.

Paul's two "optimal" results of e-13 and e-11 illustrate the dangers.

Why are most of the params and variable initial values so small?

General-purpose optimizers can't tell if small numbers are meaningful

or if they're simply noise.

In brief, choose your units so that the likely value for each variable

is of order 1 (not e-6).  Typical gradient values should also be of order 1.

You can't rely on an automatic scaling routine to arrange this.

(There are often groups of variables of similar type.  I mean typical values

for variables and gradients in any such group.)

In case Paul's comments aren't already clear, the objective values

e-13 and e-11 are probably both meaningless.  But if everything is

multiplied by e+6, there's hope for finding a reasonably accurate (local)

minimizer.

Michael

--
You received this message because you are subscribed to the Google Groups "AMPL Modeling Language" group.
To post to this group, send email to am...@googlegroups.com.
To unsubscribe from this group, send email to ampl+uns...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/ampl?hl=en.

[David Veerasingam]

Thanks Michael, my thoughts exactly -- that was the first thing that
jumped out at me. Here's my reformulation, with a few notes:

0. Scaling is very important -- it is often a good idea to scale your
variables and equations such that their ranges is approximately within
the [1,1000] range. Otherwise, you'll hit all kinds of numerical
conditioning problems (your Jacobians become ill-conditioned), and
solvers may have a hard time converging (or will converge to a
meaningless answer).
1. Variable scaling: All lengths
(MpW,MpL,MnW,MnL,minWidth,minLength,maxWidth,maxLength,LdNmos,LdPmos)
have been multiplied by 1e6 -- they are now in "mega" units. All other
variables retain their original units.
2. Equation scaling: equations with small LHS and RHS have been
premultiplied on both sides. (see constraints c1, c2, EqualCurrent)
3. Bound constraints have been rewritten as variable bounds, ie. var
x, >= U, <= L. Some solvers treat simple bounds specially, so this
formulation will aid those solvers.
4. Dummy variables (dummy1,dummy2) were introduced to remove the
divisions (quotients). Divisions are bad because they give rise to
more nonlinear Jacobian elements, and there is sometimes a possibility
of divisions by 0 during iterations.
5. This looks like a bilinear program. You can get pretty tight convex
over and underestimators for these to guide the solvers. The tightest
possible relaxation for bilinear terms (if you know the variable
bounds) are McCormick relaxations.

# ------------------------------------------------------
# Constants
param LdNmos := 4e-02; # was 4e-08
param LdPmos := 5e-02; # was 5e-08
param minWidth := 4e-01; # was 4e-07
param minLength := 3.5e-01; # was 3.5e-07
param maxWidth := 2e2; # was 0.0002
param maxLength := 2e1; # was 2e-05

param out := 1.65;
param MpCurFac := 1.144981598e-05;
param MpVth := 0.6798;
param MnCurFac := 6.455286683e-05;
param MnVth := 0.5591;
param Vnbias := 1;
param Vdd := 3.3;

# Variables
var MpW, >= minWidth, <= maxWidth := 1; # was 1e-06
var MpL, >= minLength, <= maxLength := 1; # was 1e-06
var MnW, >= minWidth, <= maxWidth := 1; # was 1e-06
var MnL, >= minLength, <= maxLength := 1; # was 1e-06
var Vpbias := 1.353;
var MnId;
var MpId;
var dummy1;
var dummy2;
# Objective
minimize area: MnW*MnL+MpW*MpL;
subject to
# Constraints
c1: 1e6 * MnId = 1e6 * MnCurFac*dummy1*(Vnbias-MnVth)^2;
c2: 1e6 * MpId = 1e6 * MpCurFac*dummy2*(Vdd-Vpbias-MpVth)^2;
c3: dummy1*(MnL-2*LdNmos) = MnW;
c4: dummy2*(MpL-2*LdPmos) = MpW;

SatP: out <= Vpbias + MpVth;

SatP2: Vdd - Vpbias >= MpVth;
EqualCurrent: 1e6 * MnId = 1e6 * MpId;
# Solve
option solver couenne;
solve;
printf "startWriteVars\n";
display area,MpW, MpL, MnW, MnL, Vpbias, MnId, MpId, dummy1, dummy2;

printf "endWriteVars\n";
display area.result;
display area.result_num;
display area.message;
display area.exitcode;

# ------------------------------------------------------

I invoked the Couenne solver (a global solver) and I managed to get a
global solution. (Incidentally, Couenne failed on the unscaled problem
-- I guess the convex envelopes must have been near infeasible because
the numbers were so close to the numerical tolerance). Here's the
output:

------------------------------------------------------
couenne: ANALYSIS TEST: Problem size before reformulation: 9 variables
(0 integer), 5 constraints.
Problem size after reformulation: 19 variables (0 integer), 1
constraints.

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

Cbc0012I Integer solution of 0.28 found by Init Rounding NLP after 0
iterations and 0 nodes (-0.00 seconds)
Cbc0001I Search completed - best objective 0.2800000037782208, took 0
iterations and 0 nodes (0.00 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost

"Finished"

couenne: Optimal
startWriteVars
area = 0.28
MpW = 0.4
MpL = 0.35
MnW = 0.4
MnL = 0.35
Vpbias = 1.61284
MnId = 1.85905e-05
MpId = 1.85905e-05
dummy1 = 1.48148
dummy2 = 1.6

endWriteVars
area.result = solved

area.result_num = 3

area.message = '\
couenne: Optimal'

area.exitcode = 0
# ------------------------------------------------------

Here are results from some local NLP solvers:
# ------------------------------------------------------
Ipopt 3.8.2:
area = 0.28
MpW = 0.4
MpL = 0.35
MnW = 0.4
MnL = 0.35
Vpbias = 1.61284
MnId = 1.85905e-05
MpId = 1.85905e-05
dummy1 = 1.48148
dummy2 = 1.6
# ------------------------------------------------------
SNOPT 7.2-8
area = 0.28
MpW = 0.4
MpL = 0.35
MnW = 0.4
MnL = 0.35
Vpbias = 1.61284
MnId = 1.85905e-05
MpId = 1.85905e-05
dummy1 = 1.48148
dummy2 = 1.6
# ------------------------------------------------------
CONOPT 3.14G
area = 0.28
MpW = 0.4
MpL = 0.35
MnW = 0.4
MnL = 0.35
Vpbias = 1.61284
MnId = 1.85905e-05
MpId = 1.85905e-05
dummy1 = 1.48148
dummy2 = 1.6


So it looks like all the solvers managed to get to the global solution
(area = 0.28 megaunits = 2.8e-7 in original units) from the given
initial guess. Lucky...

But the only way to be sure is to check it with global solvers like
Couenne or LaGO. (both free and available as AMPL solvers).


David

> > ampl+uns...@googlegroups.com <ampl%2Bunsu...@googlegroups.com>.

[David Veerasingam]

Sorry I got sidetracked, but to answer your original question....

What I usually do is to match the optimizer's output with a regexp:

if match(solve_message,'[Oo]ptimal') == 0 then {
# do something
}

The string to match varies from solver to solver. The reason I take
this approach is because some solvers don't necessary return the
correct exit codes etc. but they always return a solve_message which I
can parse.

David

[Michael Saunders]

OK David, that's a big improvement *except* for the occurrences of

1e6*MnId and 1e6*MpId (involving variables).  As you see, the final

value of those variables is 1e-5, which may or may not be noise.

Redefine the units to be millions and the optimizers will be delighted.

(You might as well redefine MnCurFac and MpCurFac too.)

Certainly it's good to avoid the divisions.

In haste,

Michael

To unsubscribe from this group, send email to ampl+uns...@googlegroups.com.

[Pietro Belotti]

Couenne indeed is sensitive to large constants -- it often happens
that with constraints and objectives containing large (>=1e6)
constants the lower bounding procedure produces a numerically instable
LP. Whether this is why Couenne can't solve the unscaled problem is a
good question, though -- but certainly it helps removing large
scalars. As for some of the constraints:

> c1: 1e6 * MnId = 1e6 * MnCurFac*dummy1*(Vnbias-MnVth)^2;
> c2: 1e6 * MpId = 1e6 * MpCurFac*dummy2*(Vdd-Vpbias-MpVth)^2;

> EqualCurrent: 1e6 * MnId = 1e6 * MpId;

The way Couenne manipulates these constraints is by creating (say for
c1) new variables y, z, and t such that

y = 1e6*MnId
z = 1e6*MnCurFac*t ~ 6.455286683*t
t = dummy1*(Vnbias-MnVth)^2

plus other variables to handle the expression defining t. The
resulting LP relaxation will contain linear constraints

y = 1e6*MnId
z = 6.455286683*t

so from Couenne's point of view, adding those 1e6 gets rid of a small
constant while introducing another (big) one. This might make a
difference for c1 and c2, but for EqualCurrent I'd doubt that. In fact
I'd simply write EqualCurrent as

EqualCurrent: MnId = MpId;

as Couenne will copy it directly to the LP relaxation.

Pietro

On Nov 16, 8:01 pm, Michael Saunders <saund...@stanford.edu> wrote:
> OK David, that's a big improvement *except* for the occurrences of
> 1e6*MnId and 1e6*MpId (involving variables).  As you see, the final
> value of those variables is 1e-5, which may or may not be noise.
> Redefine the units to be millions and the optimizers will be delighted.
> (You might as well redefine MnCurFac and MpCurFac too.)
>
> Certainly it's good to avoid the divisions.
>
> In haste,
> Michael
>

> >         For more information visithttp://projects.coin-or.org/Ipopt

> ...
>
> read more »

[brian]

Hi Michael,

This is an electrical engineering problem for integrated circuit (IC)
design. ICs are the little black boxes that make up electronics. The
numbers are so small because these are the sizes of the devices that
are used in ICs. From my master's research I knew there were scaling
issues. Scaling is on my todo list of things to fix for my PhD
research.


Hi David,

Thank you for reworking this problem!!! I wasn't sure the best way to
scale the problem. This has saved me a lot of time.


Hi Pietro,

Thank you for sharing your knowledge of Couenne.

Brian