Solving Nonlinear Programs with Piecewise‑Linear Approximation in Gurobi

Constrained optimization seeks real‑valued decision variables that satisfy a set of equality and inequality constraints while optimizing a linear or nonlinear objective. When the objective or constraints are nonlinear (NLP), solving the model directly can be difficult, but linear programming (LP) solvers such as Gurobi are fast and numerically stable.

One practical approach is to replace each nonlinear term with a piecewise‑linear (PWL) function defined by a set of breakpoints. The resulting model becomes a separable linear program that can be handled by any LP/MIP solver. The article first introduces separable functions, then presents a classic portfolio‑selection problem (Bradley et al., 1977) that is originally non‑separable because of a quadratic term.

“If the objective function is nonlinear and/or the feasible region is defined by nonlinear constraints, the problem is a nonlinear programming problem (NLP).” – Bradley et al., 1977, p. 410

By introducing an auxiliary variable and rewriting the quadratic term as a convex combination of breakpoint values, the problem is converted into a separable form. Each nonlinear term f_i(x_i) is approximated by a PWL function built from line segments whose endpoints lie on the original curve. The convex combination coefficients θ_{i,k} satisfy two conditions: they are non‑negative, sum to one, and at most two adjacent coefficients can be non‑zero. This adjacency requirement is naturally expressed with SOS‑type 2 constraints, which Gurobi can detect and exploit.

Special Ordered Set (SOS) Type 2

An SOS‑type 2 set contains variables that may be non‑zero only in pairs of adjacent indices. In Gurobi the constraint is added with model.addSOS(GRB.SOS_TYPE2, [θ_{i,0}, …, θ_{i,r‑1}]). Using SOS‑type 2 eliminates the need for extra binary variables to enforce adjacency.

Python Implementation

The implementation proceeds as follows:

Select lower and upper bounds for each decision variable.

Choose a number r of uniformly spaced breakpoints a_i = np.linspace(lb, ub, r).

Evaluate each nonlinear function at the breakpoints to obtain f_i_hat_r.

Create a Gurobi model and add continuous decision variables x_i and auxiliary variables f_i_hat.

Add variables θ_{i,k} for the convex combination (bounds 0‑1).

Link each x_i to its breakpoints with x_i == Σ θ_{i,k}·a_k and each f_i_hat to the function values with f_i_hat == Σ θ_{i,k}·f_i_hat_r[k].

Enforce the convex‑combination sum Σ θ_{i,k} == 1 for every i.

Add the original linear constraints (e.g., x0 + x1 ≤ 5, x0 + x1 - x2 == 0).

Insert SOS‑type 2 constraints for each θ vector.

Set the objective to maximize the sum of the approximated functions f0_hat + f1_hat + f2_hat and call model.optimize().

# Imports
import gurobipy as grb
import numpy as np

# Nonlinear functions
f0 = lambda x0: 20.0*x0 - 2.0*x0**2
f1 = lambda x1: 16.0*x1 - x1**2
f2 = lambda x2: -x2**2

lb, ub = 0.0, 5.0
r = 4  # number of breakpoints
a_i = np.linspace(start=lb, stop=ub, num=r)

f0_hat_r = f0(a_i)
f1_hat_r = f1(a_i)
f2_hat_r = f2(a_i)

model = grb.Model()
# decision variables
x0 = model.addVar(lb=0.0, ub=5.0, name='x0', vtype=grb.GRB.CONTINUOUS)
x1 = model.addVar(lb=0.0, ub=5.0, name='x1', vtype=grb.GRB.CONTINUOUS)
x2 = model.addVar(lb=0.0, ub=5.0, name='x2', vtype=grb.GRB.CONTINUOUS)

# auxiliary variables for function values
f0_hat = model.addVar(lb=0.0, ub=50.0, name='f0_hat', vtype=grb.GRB.CONTINUOUS)
f1_hat = model.addVar(lb=0.0, ub=55.0, name='f1_hat', vtype=grb.GRB.CONTINUOUS)
f2_hat = model.addVar(lb=-25.0, ub=0.0, name='f2_hat', vtype=grb.GRB.CONTINUOUS)

# convex‑combination variables
theta = model.addVars(3, r, lb=0.0, ub=1.0, name='theta', vtype=grb.GRB.CONTINUOUS)

# link x to breakpoints
model.addConstr(x0 == grb.quicksum(theta[0,i]*a_i[i] for i in range(r)))
model.addConstr(x1 == grb.quicksum(theta[1,i]*a_i[i] for i in range(r)))
model.addConstr(x2 == grb.quicksum(theta[2,i]*a_i[i] for i in range(r)))

# link function approximations
model.addConstr(f0_hat == grb.quicksum(theta[0,i]*f0_hat_r[i] for i in range(r)))
model.addConstr(f1_hat == grb.quicksum(theta[1,i]*f1_hat_r[i] for i in range(r)))
model.addConstr(f2_hat == grb.quicksum(theta[2,i]*f2_hat_r[i] for i in range(r)))

# convex‑combination sums
for i in range(3):
    model.addConstr(grb.quicksum(theta[i,j] for j in range(r)) == 1)

# original linear constraints
model.addConstr(x0 + x1 <= 5.0)
model.addConstr(x0 + x1 - x2 == 0.0)

# SOS‑type 2 constraints to enforce adjacency
model.addSOS(grb.GRB.SOS_TYPE2, [theta[0,i] for i in range(r)])
model.addSOS(grb.GRB.SOS_TYPE2, [theta[1,i] for i in range(r)])
model.addSOS(grb.GRB.SOS_TYPE2, [theta[2,i] for i in range(r)])

# objective
model.setObjective(f0_hat + f1_hat + f2_hat, grb.GRB.MAXIMIZE)
model.optimize()

Running the model with four breakpoints yields an optimal objective of 45, with variable values x0≈1.67, x1≈3.33, x2=5. Increasing the breakpoint count to 16 improves the objective to 46.33, demonstrating the trade‑off between approximation accuracy and model size.

Conclusion

Piecewise‑linear approximation combined with SOS‑type 2 constraints provides a systematic way to solve separable nonlinear programs using fast LP/MIP solvers. The method works for convex, concave, and mixed‑curvature functions, and the Python‑Gurobi implementation shown can be adapted to a wide range of practical optimization problems without requiring a dedicated nonlinear solver.