Nonlinear transportation model#
Description: Nonlinear transportation problem with Amplpy nltransd.mod, nltrans.dat, and nltrans.run
Tags: ampl-book, nonlinear, ipopt
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
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 nodesDEST: final nodes
Parameters:
supply {ORIG}: available units at originsdemand {DEST}: required units at destinationslimit {ORIG,DEST}: maximum capacity on routes between two nodesrate {ORIG,DEST}: base shipment costs per unit
Variables:
Trans {ORIG,DEST}: units to be shipped
Objective: minimize total cost (nonlinear)
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
import pandas as pd
import numpy as np
# ORIG supply data
supply_data = {"ORIG": ["GARY", "CLEV", "PITT"], "supply": [1400, 2600, 2900]}
supply_df = pd.DataFrame(supply_data).set_index("ORIG")
# DEST demand data
demand_data = {
"DEST": ["FRA", "DET", "LAN", "WIN", "STL", "FRE", "LAF"],
"demand": [900, 1200, 600, 400, 1700, 1100, 1000],
}
demand_df = pd.DataFrame(demand_data).set_index("DEST")
rate = np.array(
[
# 'FRA', 'DET', 'LAN', 'WIN', 'STL', 'FRE', 'LAF'
[39, 14, 11, 14, 16, 82, 8], # GARY
[27, 9, 12, 9, 26, 95, 17], # CLEV
[24, 14, 17, 13, 28, 99, 20], # PITT
]
)
limit = np.array(
[
# 'FRA', 'DET', 'LAN', 'WIN', 'STL', 'FRE', 'LAF'
[500, 1000, 1000, 1000, 800, 500, 1000], # GARY
[500, 800, 800, 800, 500, 500, 1000], # CLEV
[800, 600, 600, 600, 500, 500, 900], # PITT
]
)
ampl.read("nltrans.mod")
ampl.set_data(supply_df, "ORIG")
ampl.set_data(demand_df, "DEST")
ampl.param["rate"] = rate
ampl.param["limit"] = limit
# ipopt output is very verbose, we can reduce it by setting the option ipopt_options='outlev=0, or simply verbose=False
ampl.solve(solver="ipopt", ipopt_options="outlev=0")
ampl.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
;
From the output of the solver we can see
“Ipopt 3.12.13: Maximum Number of Iterations Exceeded.”
Since Ipopt is a local search solver, it takes several iterations to improve the current solution. This is typical in non linear problems, and this particular example is highly nonlinear. Ipopt finished the solving process after reaching the maximum iterations limit. The “solve status” can be checked through the solve_result attribute:
print(f"Solve status: {ampl.solve_result}")
assert ampl.solve_result == "limit", ampl.solve_result
Solve status: limit
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 the option send_statuses to 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.option["send_statuses"] = 0
ampl.eval("let {i in ORIG, j in DEST} Trans[i,j] := limit[i,j]/2;")
ampl.solve(solver="ipopt", verbose=False)
ampl.display("Total_Cost")
ampl.display("{i in ORIG, j in DEST} Trans[i,j]/limit[i,j]")
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
;
print(f"Solve status: {ampl.solve_result}")
assert ampl.solve_result == "solved", ampl.solve_result
Solve status: solved
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:
%%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.reset()
ampl.read("nltransc.mod")
ampl.set_data(supply_df, "ORIG")
ampl.set_data(demand_df, "DEST")
ampl.param["rate"] = rate
ampl.param["limit"] = limit
ampl.solve(solver="ipopt", ipopt_options="outlev=0")
ampl.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>
Total_Cost = 0
Trans [*,*] (tr)
: CLEV GARY PITT :=
DET 0 0 0
FRA 0 0 0
FRE 0 0 0
LAF 0 0 0
LAN 0 0 0
STL 0 0 0
WIN 0 0 0
;
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.reset()
ampl.read("nltransd.mod")
ampl.set_data(supply_df, "ORIG")
ampl.set_data(demand_df, "DEST")
ampl.param["rate"] = rate
ampl.param["limit"] = limit
ampl.solve(solver="ipopt", verbose=False)
ampl.display("Total_Cost", "Trans")
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.reset()
ampl.read("nltranse.mod")
ampl.set_data(supply_df, "ORIG")
ampl.set_data(demand_df, "DEST")
ampl.param["rate"] = rate
ampl.param["limit"] = limit
num_iterations = 10
for it in range(1, num_iterations + 1):
ampl.param["alpha"] = 0.1 * it
ampl.solve(solver="ipopt", verbose=False)
print(f"-- Iteration {it}")
ampl.display("alpha", "Total_Cost")
-- Iteration 1
alpha = 0.1
Total_Cost = 354287
-- Iteration 2
alpha = 0.2
Total_Cost = 354278
-- Iteration 3
alpha = 0.3
Total_Cost = 354277
-- Iteration 4
alpha = 0.4
Total_Cost = 354279
-- Iteration 5
alpha = 0.5
Total_Cost = 354278
-- Iteration 6
alpha = 0.6
Total_Cost = 354278
-- Iteration 7
alpha = 0.7
Total_Cost = 354277
-- Iteration 8
alpha = 0.8
Total_Cost = 354277
-- Iteration 9
alpha = 0.9
Total_Cost = 354278
-- Iteration 10
alpha = 1
Total_Cost = 354278