Linear Programming with GTSAM

This notebook solves a two-dimensional LP with one equality constraint and active inequality bounds. Plotly contours and constraint traces show how the active bound selects the solution on the feasible line segment.

import numpy as np
import gtsam
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from IPython.display import display
from gtsam.symbol_shorthand import X

np.set_printoptions(precision=4, suppress=True)

A 2D LP with an active bound¶

We minimize -x, so the optimizer pushes right along the equality line x + y = 1 until the upper bound x <= 0.35 becomes active.

key = X(0)
model = gtsam.noiseModel.Unit.Create(1)

problem = gtsam.LpProblem()
problem.addCost(gtsam.JacobianFactor(key, np.array([[-1.0, 0.0]]), np.zeros(1), model))
problem.addConstraint(
    gtsam.LinearConstraint.Equal(
        gtsam.JacobianFactor(key, np.array([[1.0, 1.0]]), np.array([1.0]), model)
    )
)
problem.addConstraint(
    gtsam.LinearConstraint.LessEqual(
        gtsam.JacobianFactor(key, np.array([[1.0, 0.0]]), np.array([0.35]), model)
    )
)
problem.addConstraint(
    gtsam.LinearConstraint.GreaterEqual(
        gtsam.JacobianFactor(
            key, np.eye(2), np.zeros(2), gtsam.noiseModel.Unit.Create(2)
        )
    )
)
problem.addConstraint(
    gtsam.LinearConstraint.GreaterEqual(
        gtsam.JacobianFactor(key, np.array([[0.0, 1.0]]), np.array([0.65]), model)
    )
)

initial = gtsam.Values()
initial.insert(key, np.array([0.2, 0.8]))
result = problem.optimize(initial)
solution = result.atVector(key)
cost, equality_violation, inequality_violation = problem.evaluate(result)

print(f"solution = {solution}")
print(f"objective = {cost:.6f}")
print(
    f"violations: equality={equality_violation:.2e}, inequality={inequality_violation:.2e}"
)
print(f"active x upper bound residual = {solution[0] - 0.35:+.2e}")
print(f"active y lower bound residual = {0.65 - solution[1]:+.2e}")

solution = [0.35 0.65]
objective = 0.000000
violations: equality=0.00e+00, inequality=0.00e+00
active x upper bound residual = +0.00e+00
active y lower bound residual = +0.00e+00

grid_x = np.linspace(-0.1, 1.1, 160)
grid_y = np.linspace(-0.1, 1.1, 160)
X_grid, Y_grid = np.meshgrid(grid_x, grid_y)
objective = -X_grid
line_x = np.linspace(-0.05, 1.05, 120)
line_y = 1.0 - line_x
feasible_x = np.linspace(0.0, 0.35, 60)
feasible_y = 1.0 - feasible_x

fig = go.Figure()
fig.add_trace(go.Contour(x=grid_x, y=grid_y, z=objective, colorscale='Viridis', opacity=0.65, contours=dict(showlabels=True), name='linear objective'))
fig.add_trace(go.Scatter(x=line_x, y=line_y, mode='lines', line=dict(color='black', dash='dash'), name='x + y = 1'))
fig.add_trace(go.Scatter(x=feasible_x, y=feasible_y, mode='lines', line=dict(color='seagreen', width=7), name='feasible equality segment'))
fig.add_trace(go.Scatter(x=[0.35, 0.35], y=[-0.1, 1.1], mode='lines', line=dict(color='crimson', width=3), name='x <= 0.35'))
fig.add_trace(go.Scatter(x=[0.0, 1.1], y=[0.0, 0.0], mode='lines', line=dict(color='royalblue', width=2), name='y >= 0'))
fig.add_trace(go.Scatter(x=[-0.1, 1.1], y=[0.65, 0.65], mode='lines', line=dict(color='crimson', width=3, dash='dot'), name='y >= 0.65'))
fig.add_trace(go.Scatter(x=[0.0, 0.0], y=[0.0, 1.1], mode='lines', line=dict(color='royalblue', width=2), name='x >= 0'))
fig.add_trace(go.Scatter(x=[solution[0]], y=[solution[1]], mode='markers+text', marker=dict(size=15, color='white', line=dict(color='black', width=2)), text=['solution'], textposition='top center', name='solution'))
fig.update_layout(template='plotly_white', width=800, height=620, title='LP geometry: two active inequalities on the equality line', xaxis_title='x', yaxis_title='y', yaxis_scaleanchor='x', legend=dict(x=0.02, y=0.98))
display(fig)