Nonlinear transportation model#

nltrans_lecture.ipynb Open In Colab Kaggle Gradient Open In SageMaker Studio Lab Hits

Description: book example autogenerated using nltransd.mod, nltrans.dat, and nltrans.run

Tags: ampl-only, ampl-book, nonlinear

Notebook author: Marcos Dominguez Velad <marcos@ampl.com>

Model author: N/A

# Install dependencies
%pip install -q amplpy

# Google Colab & Kaggle integration
from amplpy import AMPL, ampl_notebook

ampl = ampl_notebook(
    modules=["coin"],  # modules to install
    license_uuid="default",  # license to use
)  # instantiate AMPL object and register magics

Nonlinear transportation model#

This is a variation of the linear transportation model presented on the Chapter 3 of the AMPL book, containing a nonlinear objective. There are a set of origin nodes, and a set of destination nodes (net model).

  • Sets:

    • ORIG: origin nodes

    • DEST: final nodes

  • Parameters:

    • supply {ORIG}: available units at origins

    • demand {DEST}: required units at destinations

    • limit {ORIG,DEST}: maximum capacity on routes between two nodes

    • rate {ORIG,DEST}: base shipment costs per unit

  • Variables:

    • Trans {ORIG,DEST}: units to be shipped

  • Objective: minimize total cost (nonlinear)

\[\begin{split}\sum \limits_{\substack{i \in ORIG \\ j \in DEST}} rate[i,j] \cdot \frac{Trans[i,j]}{1 - \frac{Trans[i,j]}{limit[i,j]}}\end{split}\]

The bigger the Trans[i,j] value, the closer to limit[i,j] (upper bound) so denominator tends to \(1-1=0\) implying high costs.

  • Constraints:

    • Supply {ORIG}: node ships units equal to supply capacity:

    \[\sum \limits_{j \in DEST} Trans[i,j] = supply[i]\]

    • Demand {DEST}: node gets units equal to demand:

    \[\sum \limits_{i \in ORIG} Trans[i,j] = demand[j]\]
%%writefile nltrans.mod
set ORIG;   # origins
set DEST;   # destinations

param supply {ORIG} >= 0;   # amounts available at origins
param demand {DEST} >= 0;   # amounts required at destinations

   check: sum {i in ORIG} supply[i] = sum {j in DEST} demand[j];

param rate {ORIG,DEST} >= 0;   # base shipment costs per unit
param limit {ORIG,DEST} > 0;   # limit on units shipped

var Trans {i in ORIG, j in DEST} >= 0;    # actual units to be shipped

minimize Total_Cost:
   sum {i in ORIG, j in DEST}
      rate[i,j] * Trans[i,j] / (1 - Trans[i,j]/limit[i,j]);

subject to Supply {i in ORIG}:  
   sum {j in DEST} Trans[i,j] = supply[i];

subject to Demand {j in DEST}:  
   sum {i in ORIG} Trans[i,j] = demand[j];
   

Overwriting nltrans.mod

%%writefile nltrans.dat
data;

param: ORIG:  supply :=
        GARY   1400    CLEV   2600    PITT   2900 ;

param: DEST:  demand :=
        FRA     900    DET    1200    LAN     600 
        WIN     400    STL    1700    FRE    1100 
        LAF    1000 ;

param rate :  FRA  DET  LAN  WIN  STL  FRE  LAF :=
        GARY   39   14   11   14   16   82    8
        CLEV   27    9   12    9   26   95   17
        PITT   24   14   17   13   28   99   20 ;

param limit :  FRA  DET  LAN  WIN  STL  FRE  LAF :=
        GARY   500 1000 1000 1000  800  500 1000
        CLEV   500  800  800  800  500  500 1000
        PITT   800  600  600  600  500  500  900 ;

Overwriting nltrans.dat

%%ampl_eval
model nltrans.mod;
data nltrans.dat;
option solver ipopt;
option ipopt_options 'outlev=0';
solve;
display Total_Cost, Trans;

Ipopt 3.12.13: outlev=0


******************************************************************************
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 http://projects.coin-or.org/Ipopt
******************************************************************************

 
Ipopt 3.12.13: Maximum Number of Iterations Exceeded.

suffix ipopt_zU_out OUT;
suffix ipopt_zL_out OUT;
Total_Cost = -22681400000

Trans [*,*] (tr)
:      CLEV       GARY      PITT      :=
DET   518.41    279.922    401.668
FRA   326.573   153.849    419.579
FRE   377.651   371.911    350.438
LAF   409.194   173.889    416.917
LAN   284.738    34.4686   280.794
STL   500       372.77     827.229
WIN   183.433    13.1917   203.375
;

assert ampl.solve_result in ["solved", "limit"], ampl.solve_result

The previous problem is unbounded, however, this is not what we expect for this problem. This is due to the highly nonlinear objective function.

Because the behavior of a nonlinear optimization algorithm can be sensitive to the choice of starting guess, we might hope that the solver will have greater success from a different start. To ensure that the comparison is meaningful, we first set

option send_statuses 0;

So that status information about variables that was returned by the previous solve will not be sent back to help determine a starting point for the next solve. We can use the let command to suggest a new initial value for each variable. For example, say that Trans[i,j] is half of limit[i,j].

%%ampl_eval
option send_statuses 0;
let {i in ORIG, j in DEST} Trans[i,j] := limit[i,j]/2;
solve;
display Total_Cost;
display {i in ORIG, j in DEST} Trans[i,j]/limit[i,j];

Ipopt 3.12.13: outlev=0


******************************************************************************
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 http://projects.coin-or.org/Ipopt
******************************************************************************

 
Ipopt 3.12.13: Optimal Solution Found
Total_Cost = 1212120

Trans[i,j]/limit[i,j] [*,*] (tr)
:       CLEV       GARY       PITT      :=
DET   0.732965   0.187385   0.710404
FRA   0.589986   0.162441   0.654733
FRE   0.731      0.739444   0.729557
LAF   0.490537   0          0.56607
LAN   0.367685   0          0.509753
STL   0.939382   0.95209    0.937274
WIN   0.123449   0          0.502068
;

assert ampl.solve_result == "solved", ampl.solve_result

Now the solver gives a right solution. However, we could have had trouble in case any Trans variable get closer to the limit, what we could check by calling

display {i in ORIG, j in DEST} Trans[i,j]/limit[i,j];

We could implement a more robust formulation by giving more explicit bounds to the variables to avoid numerical issues. With the following bounds, Trans will be smaller that limit.

var Trans {i in ORIG, j in DEST} >= 0, <= .9999 * limit[i,j];

Multiple local optima#

Let’s change to a concave objective function by adding a power of 0.8:

\[\begin{split}\sum \limits_{\substack{i \in ORIG \\ j \in DEST}} rate[i,j] \cdot \frac{Trans[i,j]^{0.8}}{1 - \frac{Trans[i,j]}{limit[i,j]}}\end{split}\]
%%writefile nltransc.mod
set ORIG;   # origins
set DEST;   # destinations

param supply {ORIG} >= 0;   # amounts available at origins
param demand {DEST} >= 0;   # amounts required at destinations

   check: sum {i in ORIG} supply[i] = sum {j in DEST} demand[j];

param rate {ORIG,DEST} >= 0;   # base shipment costs per unit
param limit {ORIG,DEST} > 0;   # limit on units shipped

var Trans {i in ORIG, j in DEST} >= 0, <= .9999 * limit[i,j];
    # actual units to be shipped

minimize Total_Cost:
   sum {i in ORIG, j in DEST}
      rate[i,j] * Trans[i,j]^0.8 / (1 - Trans[i,j]/limit[i,j]);

subject to Supply {i in ORIG}:  
   sum {j in DEST} Trans[i,j] = supply[i];

subject to Demand {j in DEST}:  
   sum {i in ORIG} Trans[i,j] = demand[j];

Overwriting nltransc.mod

%%ampl_eval
reset;
model nltransc.mod;
data nltrans.dat;
option solver ipopt;
option ipopt_options 'outlev=0';
solve;
display Total_Cost, Trans;

Ipopt 3.12.13: outlev=0


******************************************************************************
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 http://projects.coin-or.org/Ipopt
******************************************************************************

Error evaluating "var =" definition -1: can't evaluate pow'(0,0.8).
exit value 1
<BREAK>

assert ampl.solve_result == "failure", ampl.solve_result

We got the error

Error evaluating “var =” definition -1: can’t evaluate pow’(0,0.8).

This is due to the fact that the derivative of the exponential function at zero is infinity. We could just give an starting point away from zero, and tighten the lower bound too:

var Trans {i in ORIG, j in DEST} >= 1e-10, <= .9999 * limit[i,j], := limit[i,j]/2;

%%writefile nltransd.mod
set ORIG;   # origins
set DEST;   # destinations

param supply {ORIG} >= 0;   # amounts available at origins
param demand {DEST} >= 0;   # amounts required at destinations

   check: sum {i in ORIG} supply[i] = sum {j in DEST} demand[j];

param rate {ORIG,DEST} >= 0;   # base shipment costs per unit
param limit {ORIG,DEST} > 0;   # limit on units shipped

var Trans {i in ORIG, j in DEST} >= 1e-10, <= .9999 * limit[i,j], := limit[i,j]/2;
    # actual units to be shipped

minimize Total_Cost:
   sum {i in ORIG, j in DEST}
      rate[i,j] * Trans[i,j]^0.8 / (1 - Trans[i,j]/limit[i,j]);

subject to Supply {i in ORIG}:  
   sum {j in DEST} Trans[i,j] = supply[i];

subject to Demand {j in DEST}:  
   sum {i in ORIG} Trans[i,j] = demand[j];

Overwriting nltransd.mod

%%ampl_eval
reset;
model nltransd.mod;
data nltrans.dat;
option solver ipopt;
option ipopt_options 'outlev=0';
solve;
display Total_Cost, Trans;

Ipopt 3.12.13: outlev=0


******************************************************************************
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 http://projects.coin-or.org/Ipopt
******************************************************************************

 
Ipopt 3.12.13: Restoration Phase Failed.

suffix ipopt_zU_out OUT;
suffix ipopt_zL_out OUT;
Total_Cost = 354279

Trans [*,*] (tr)
:      CLEV         GARY         PITT      :=
DET   586.429   191.805         421.765
FRA   292.089    75.1764        532.735
FRE   365.333   370.173         364.494
LAF   488.914     0.0937081     510.992
LAN   298.74      0.000276411   301.26
STL   469.177   762.751         468.072
WIN    99.318     9.99922e-05   300.682
;

assert ampl.solve_result == "failure", ampl.solve_result

By adding a new parameter alpha that ponders the lower bound for Trans, we may solve the problem from different initial guesses, providing potential different solutions.

%%writefile nltranse.mod
set ORIG;   # origins
set DEST;   # destinations

param supply {ORIG} >= 0;   # amounts available at origins
param demand {DEST} >= 0;   # amounts required at destinations

   check: sum {i in ORIG} supply[i] = sum {j in DEST} demand[j];

param rate {ORIG,DEST} >= 0;   # base shipment costs per unit
param limit {ORIG,DEST} > 0;   # limit on units shipped

param alpha >= 0, <= 1;
var Trans {i in ORIG, j in DEST} >= 1e-10, <= .9999 * limit[i,j], := alpha*limit[i,j];
    # actual units to be shipped

minimize Total_Cost:
   sum {i in ORIG, j in DEST}
      rate[i,j] * Trans[i,j]^0.8 / (1 - Trans[i,j]/limit[i,j]);

subject to Supply {i in ORIG}:  
   sum {j in DEST} Trans[i,j] = supply[i];

subject to Demand {j in DEST}:  
   sum {i in ORIG} Trans[i,j] = demand[j];

Overwriting nltranse.mod

%%ampl_eval
reset;
set TOTAL_VALUES := 1..10;
model nltranse.mod;
data nltrans.dat;
option solver ipopt;
option ipopt_options 'outlev=0';
for {it in TOTAL_VALUES}{
    let alpha := 0.1*it;
    solve;
    printf "## Iteration %d\n", it;
    display alpha, Total_Cost;
}

Ipopt 3.12.13: outlev=0


******************************************************************************
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 http://projects.coin-or.org/Ipopt
******************************************************************************

 
Ipopt 3.12.13: Search Direction becomes Too Small.

suffix ipopt_zU_out OUT;
suffix ipopt_zL_out OUT;
## Iteration 1
alpha = 0.1
Total_Cost = 354281

Ipopt 3.12.13: outlev=0


******************************************************************************
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 http://projects.coin-or.org/Ipopt
******************************************************************************

 
Ipopt 3.12.13: Search Direction becomes Too Small.
## Iteration 2
alpha = 0.2
Total_Cost = 354277

Ipopt 3.12.13: outlev=0


******************************************************************************
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 http://projects.coin-or.org/Ipopt
******************************************************************************

 
Ipopt 3.12.13: Search Direction becomes Too Small.
## Iteration 3
alpha = 0.3
Total_Cost = 354277

Ipopt 3.12.13: outlev=0


******************************************************************************
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 http://projects.coin-or.org/Ipopt
******************************************************************************

 
Ipopt 3.12.13: Restoration Phase Failed.
## Iteration 4
alpha = 0.4
Total_Cost = 354277

Ipopt 3.12.13: outlev=0


******************************************************************************
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 http://projects.coin-or.org/Ipopt
******************************************************************************

 
Ipopt 3.12.13: Search Direction becomes Too Small.
## Iteration 5
alpha = 0.5
Total_Cost = 354278

Ipopt 3.12.13: outlev=0


******************************************************************************
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 http://projects.coin-or.org/Ipopt
******************************************************************************

 
Ipopt 3.12.13: Restoration Phase Failed.
## Iteration 6
alpha = 0.6
Total_Cost = 354278

Ipopt 3.12.13: outlev=0


******************************************************************************
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 http://projects.coin-or.org/Ipopt
******************************************************************************

 
Ipopt 3.12.13: Search Direction becomes Too Small.
## Iteration 7
alpha = 0.7
Total_Cost = 354277

Ipopt 3.12.13: outlev=0


******************************************************************************
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 http://projects.coin-or.org/Ipopt
******************************************************************************

 
Ipopt 3.12.13: Search Direction becomes Too Small.
## Iteration 8
alpha = 0.8
Total_Cost = 354278

Ipopt 3.12.13: outlev=0


******************************************************************************
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 http://projects.coin-or.org/Ipopt
******************************************************************************

 
Ipopt 3.12.13: Restoration Phase Failed.
## Iteration 9
alpha = 0.9
Total_Cost = 354278

Ipopt 3.12.13: outlev=0


******************************************************************************
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 http://projects.coin-or.org/Ipopt
******************************************************************************

 
Ipopt 3.12.13: Search Direction becomes Too Small.
## Iteration 10
alpha = 1
Total_Cost = 354277

assert ampl.solve_result == "solved", ampl.solve_result