Nonlinear transportation problem example#

Description: book example autogenerated using 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 pandas numpy

# 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{array}{r} \sum_{\begin{array}{c} i \in O R I G \\ j \in D E S T \end{array}} r a t e [i, j] \cdot \frac{T r a n s [i, j]^{0.8}}{1 - \frac{T r a n s [i, j]}{l i m i t [i, j]}} \end{array}

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_{j \in D E S T} T r a n s [i, j] = s u p p l y [i]$
- Demand {DEST}: node gets units equal to demand:
$\sum_{i \in O R I G} T r a n s [i, j] = d e m a n d [j]$

Remark: as the objective function is highly nonlinear, some bounds are fixed in order to help the solver and getting a consistent solution. The first guess for variables is away from zero in 0.5*limit[i,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}
   >= 1e-10, <= .9999 * limit[i,j], := 0.5*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 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
    ]
)

# Read the model
ampl.reset()
ampl.read("nltrans.mod")

# Load the data
ampl.set_data(supply_df, "ORIG")
ampl.set_data(demand_df, "DEST")
ampl.param["rate"] = rate
ampl.param["limit"] = limit

# ampl.display('ORIG, DEST, supply, demand, rate, limit')

# Solve the problem
# Set iterations limit and time limit for ipopt
ampl.option["ipopt_options"] = "max_iter=100 max_cpu_time=10"
ampl.solve(solver="ipopt", verbose=False)

# Show solution
Total_Cost = ampl.obj["Total_Cost"].value()
print(f"Total cost: {Total_Cost}")

Trans = ampl.var["Trans"].to_pandas()
print(Trans)

Total cost: 354280.4809582299
                Trans.val
index0 index1            
CLEV   DET     586.432409
       FRA     292.093029
       FRE     365.334291
       LAF     488.917726
       LAN     298.728286
       STL     469.182376
       WIN      99.311883
GARY   DET     191.798629
       FRA      75.168128
       FRE     370.169674
       LAF       0.086271
       LAN       0.022464
       STL     762.744843
       WIN       0.009992
PITT   DET     421.768962
       FRA     532.738843
       FRE     364.496035
       LAF     510.996003
       LAN     301.249250
       STL     468.072781
       WIN     300.678125

The solver ipopt reached maximum number of iterations while searching for the optimal solution, so the status of the problem is not ‘solved’, but ‘limit’. This can be checked through the attribute: ampl.solve_result

assert ampl.solve_result == "limit", ampl.solve_result