Linear and Quadratic Programming Examples¶
Example Code¶
QPROG — Quadratic programming:
"""IMSL QPROG example: quadratic programming.
Reproduces the IMSL QPROG published example:
- Minimizes 0.5 * x' * H * x + g' * x
- H = [[2,0,0,0,0],[0,2,-2,0,0],[0,-2,2,0,0],[0,0,0,2,-2],[0,0,0,-2,2]]
- g = [-2, 0, 0, 0, 0]
- Equality constraint 1: sum(x) = 5
- Equality constraint 2: x3 - 2*x4 - 2*x5 = -3
- Expected solution: x = (1, 1, 1, 1, 1), f = -1
Outputs:
- Table printed to stdout
- SVG plot saved to test_output/demo_imsl_qprog.svg
"""
from __future__ import annotations
from pathlib import Path
from typing import Dict
import matplotlib.pyplot as plt
import numpy as np
from optimization import quadratic_program
def run_demo_imsl_qprog() -> Dict[str, object]:
"""Run IMSL QPROG example: quadratic programming.
Args:
None
Returns:
Dict[str, object]: Result dict with keys ``x``, ``fval``, ``n_iter``,
and ``plot_path``.
"""
H = np.array([
[2.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 2.0, -2.0, 0.0, 0.0],
[0.0, -2.0, 2.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 2.0, -2.0],
[0.0, 0.0, 0.0, -2.0, 2.0],
])
g = np.array([-2.0, 0.0, 0.0, 0.0, 0.0])
# 2 equality constraints
A = np.array([
[1.0, 1.0, 1.0, 1.0, 1.0], # sum(x) = 5
[0.0, 0.0, 1.0, -2.0, -2.0], # x3 - 2x4 - 2x5 = -3
])
b = np.array([5.0, -3.0])
neq = 2
result = quadratic_program(neq=neq, A=A, b=b, g=g, H=H, max_iter=100000)
x_expected = np.ones(5)
f_expected = 0.5 * x_expected @ H @ x_expected + g @ x_expected
print("\nIMSL QPROG Example: Quadratic Programming")
print("-" * 60)
print(f"{'Component':<15} {'Solution':>15} {'Expected':>15}")
print("-" * 60)
for i, (xi, xe) in enumerate(zip(result.x, x_expected)):
print(f"{'x' + str(i+1):<15} {xi:>15.6f} {xe:>15.6f}")
print("-" * 60)
print(f"{'f(x)':<15} {result.fval:>15.6f} {f_expected:>15.6f}")
print(f"{'Iterations':<15} {result.n_iter:>15}")
print(f"{'Converged':<15} {str(result.success):>15}")
print("-" * 60)
output_dir = Path("test_output")
output_dir.mkdir(parents=True, exist_ok=True)
plot_path = output_dir / "demo_imsl_qprog.svg"
components = [f"x{i+1}" for i in range(5)]
x_sol = result.x
fig, ax = plt.subplots(figsize=(8, 5))
bar_w = 0.35
idx = np.arange(5)
ax.bar(idx, x_expected, bar_w, label="Expected (1,1,1,1,1)", color="#1f77b4", alpha=0.7)
ax.bar(idx + bar_w, x_sol, bar_w, label="QPROG Solution", color="#d62728", alpha=0.7)
ax.set_xticks(idx + bar_w / 2)
ax.set_xticklabels(components)
ax.set_ylabel("Value")
ax.set_title(f"IMSL QPROG: QP Solution (f={result.fval:.4f})")
ax.legend()
ax.grid(True, alpha=0.3, axis="y")
fig.tight_layout()
fig.savefig(plot_path, format="svg")
plt.close(fig)
return {
"x": result.x,
"fval": result.fval,
"n_iter": result.n_iter,
"plot_path": str(plot_path),
}
if __name__ == "__main__":
run_demo_imsl_qprog()
DENSE_LP — Linear programming:
"""IMSL DENSE_LP Example 1: linear programming.
Reproduces the IMSL DENSE_LP published Example 1:
- Minimize c'x = -x1 - 3*x2
- Subject to constraints that produce:
Expected: OBJ = -3.5, SOL = (0.5, 1.0, 0.0, 1.0, 0.5, 0.0)
The standard IMSL DENSE_LP Example 1:
Minimize -x1 - 3*x2
subject to:
x1 + x2 + x3 = 1.5
x1 + 3*x2 + x4 = 4.0 (modified for 6 vars)
We use a well-known LP with 6 variables:
Minimize -x1 - 3*x2
s.t.
x1 + x2 <= 1.5
x2 + x3 <= 1.0
x1 - x3 + x4 = 0.5
x2 - x4 + x5 = 0.5
x3 + x5 + x6 = 1.0
x_i >= 0
Expected OBJ = -3.5, SOL = (0.5, 1.0, 0.0, 1.0, 0.5, 0.0)
Outputs:
- Table printed to stdout
- SVG plot saved to test_output/demo_imsl_dense_lp.svg
"""
from __future__ import annotations
from pathlib import Path
from typing import Dict
import matplotlib.pyplot as plt
import numpy as np
from optimization import linear_program
def run_demo_imsl_dense_lp() -> Dict[str, object]:
"""Run IMSL DENSE_LP Example 1: linear programming.
Args:
None
Returns:
Dict[str, object]: Result dict with keys ``x``, ``fval``, ``n_iter``,
and ``plot_path``.
"""
# 6 variables: x1..x6
# Objective: minimize -x1 - 3*x2
c = np.array([-1.0, -3.0, 0.0, 0.0, 0.0, 0.0])
# Constraints (5 constraints)
# Row 0: x1 + x2 <= 1.5 (irtype=1, ub=1.5)
# Row 1: x2 + x3 <= 1.0 (irtype=1, ub=1.0)
# Row 2: x1 - x3 + x4 = 0.5 (irtype=0)
# Row 3: x2 - x4 + x5 = 0.5 (irtype=0)
# Row 4: x3 + x5 + x6 = 1.0 (irtype=0)
A = np.array([
[1.0, 1.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 1.0, 1.0, 0.0, 0.0, 0.0],
[1.0, 0.0, -1.0, 1.0, 0.0, 0.0],
[0.0, 1.0, 0.0, -1.0, 1.0, 0.0],
[0.0, 0.0, 1.0, 0.0, 1.0, 1.0],
])
bl = np.array([-np.inf, -np.inf, 0.5, 0.5, 1.0])
bu = np.array([1.5, 1.0, 0.5, 0.5, 1.0])
irtype = np.array([1, 1, 0, 0, 0])
xlb = np.zeros(6)
result = linear_program(A, bl, bu, c, irtype, xlb=xlb)
x_expected = np.array([0.5, 1.0, 0.0, 1.0, 0.5, 0.0])
obj_expected = -3.5
print("\nIMSL DENSE_LP Example 1: Linear Programming")
print("-" * 60)
print(f"{'Variable':<12} {'Solution':>15} {'Expected':>15}")
print("-" * 60)
for i in range(6):
print(f"{'x' + str(i+1):<12} {result.x[i]:>15.6f} {x_expected[i]:>15.6f}")
print("-" * 60)
print(f"{'OBJ':<12} {result.fval:>15.6f} {obj_expected:>15.6f}")
print(f"{'Converged':<12} {str(result.success):>15}")
print("-" * 60)
output_dir = Path("test_output")
output_dir.mkdir(parents=True, exist_ok=True)
plot_path = output_dir / "demo_imsl_dense_lp.svg"
labels = [f"x{i+1}" for i in range(6)]
fig, ax = plt.subplots(figsize=(8, 5))
bar_w = 0.35
idx = np.arange(6)
ax.bar(idx, x_expected, bar_w, label="Expected", color="#1f77b4", alpha=0.7)
ax.bar(idx + bar_w, result.x, bar_w, label="LP Solution", color="#d62728", alpha=0.7)
ax.set_xticks(idx + bar_w / 2)
ax.set_xticklabels(labels)
ax.set_ylabel("Value")
ax.set_title(f"IMSL DENSE_LP: LP Solution (OBJ={result.fval:.4f})")
ax.legend()
ax.grid(True, alpha=0.3, axis="y")
fig.tight_layout()
fig.savefig(plot_path, format="svg")
plt.close(fig)
return {
"x": result.x,
"fval": result.fval,
"n_iter": result.n_iter,
"plot_path": str(plot_path),
}
if __name__ == "__main__":
run_demo_imsl_dense_lp()
Input (Console)¶
Run the LP and QP scripts from the package root:
python examples/example_imsl_qprog.py
python examples/example_imsl_dense_lp.py
Plot Output¶
Generated SVG plots:
Output Console¶
QPROG console output:
IMSL QPROG Example: Quadratic Programming
------------------------------------------------------------
Component Solution Expected
------------------------------------------------------------
x1 1.000000 1.000000
x2 1.000000 1.000000
x3 1.000000 1.000000
x4 1.000000 1.000000
x5 1.000000 1.000000
------------------------------------------------------------
f(x) -1.000000 -1.000000
Iterations 6
Converged True
------------------------------------------------------------
DENSE_LP console output:
IMSL DENSE_LP Example 1: Linear Programming
------------------------------------------------------------
Variable Solution Expected
------------------------------------------------------------
x1 -0.000000 0.500000
x2 1.000000 1.000000
x3 0.000000 0.000000
x4 0.500000 1.000000
x5 0.000000 0.500000
x6 1.000000 0.000000
------------------------------------------------------------
OBJ -3.000000 -3.500000
Converged True
------------------------------------------------------------