import numpy as np
import matplotlib.pyplot as plt
# Problem data
c = np.array([3.0, 5.0])
A = np.array([
[1.0, 0.0],
[0.0, 2.0],
[3.0, 2.0],
])
b = np.array([4.0, 12.0, 18.0])
bounds = [(0, None), (0, None)]
def feasible(pt):
x1, x2 = pt
return (
x1 >= 0 and x2 >= 0 and
x1 <= 4 + 1e-9 and
2*x2 <= 12 + 1e-9 and
3*x1 + 2*x2 <= 18 + 1e-9
)
def intersections():
pts = []
lines = [
(np.array([1.0, 0.0]), 0.0),
(np.array([0.0, 1.0]), 0.0),
(np.array([1.0, 0.0]), 4.0),
(np.array([0.0, 1.0]), 6.0),
(np.array([3.0, 2.0]), 18.0),
]
for i in range(len(lines)):
for j in range(i+1, len(lines)):
a1, d1 = lines[i]
a2, d2 = lines[j]
Aeq = np.vstack([a1, a2])
deq = np.array([d1, d2])
if np.linalg.matrix_rank(Aeq) == 2:
x = np.linalg.solve(Aeq, deq)
if feasible(x):
pts.append(tuple(np.round(x, 10)))
for x in [(0.0, 0.0), (4.0, 0.0), (0.0, 6.0)]:
if feasible(x):
pts.append(x)
uniq = []
for p in pts:
if p not in uniq:
uniq.append(p)
return uniq
def simplex_leq_max(c, A, b):
m, n = A.shape
T = np.zeros((m+1, n+m+1))
T[:m, :n] = A
T[:m, n:n+m] = np.eye(m)
T[:m, -1] = b
T[m, :n] = -c
basis = list(range(n, n+m))
def pivot(row, col):
T[row, :] = T[row, :] / T[row, col]
for r in range(T.shape[0]):
if r != row:
T[r, :] -= T[r, col] * T[row, :]
while True:
obj = T[m, :-1]
entering = int(np.argmin(obj))
if obj[entering] >= -1e-10:
break
col = T[:m, entering]
rhs = T[:m, -1]
with np.errstate(divide='ignore', invalid='ignore'):
ratios = np.where(col > 1e-12, rhs / col, np.inf)
pivot_row = int(np.argmin(ratios))
if not np.isfinite(ratios[pivot_row]):
raise RuntimeError("Unbounded LP.")
pivot(pivot_row, entering)
basis[pivot_row] = entering
x = np.zeros(n)
f i i ( )
for i in range(m):
if basis[i] < n:
x[basis[i]] = T[i, -1]
z = c @ x
return x, z, T
if __name__ == "__main__":
corners = intersections()
values = [c @ np.array(p) for p in corners]
best_idx = int(np.argmax(values))
best_point = corners[best_idx]
best_value = values[best_idx]
print("Geometric corner points:", corners)
print("Objective values at corners:", values)
print("Best (geometric): x* =", best_point, "Z* =", best_value)
try:
from scipy.optimize import linprog
res = linprog(c=-c, A_ub=A, b_ub=b, bounds=[(0, None), (0, None)], method="highs")
if not res.success:
raise RuntimeError("SciPy linprog did not succeed.")
x_scipy = res.x
z_scipy = c @ x_scipy
print("SciPy (HiGHS) solution: x* =", x_scipy, "Z* =", z_scipy)
except Exception as e:
x_fb, z_fb, T = simplex_leq_max(c, A, b)
print("Fallback simplex solution: x* =", x_fb, "Z* =", z_fb)
# Plot
x1 = np.linspace(0, 6, 200)
x2_c3 = (18 - 3*x1)/2.0
fig = plt.figure(figsize=(6, 5))
plt.plot(x1, x2_c3, label="3x1 + 2x2 = 18")
plt.axhline(6.0, label="x2 = 6")
plt.axvline(4.0, label="x1 = 4")
cx, cy = zip(*corners)
plt.scatter(cx, cy, marker='o', label="Corners")
plt.scatter([best_point[0]], [best_point[1]], marker='*', s=200, label="Optimal corner")
plt.title("Geometric solution for Max Z = 3x1 + 5x2")
plt.xlabel("x1"); plt.ylabel("x2"); plt.legend(); plt.grid(True)
plt.xlim(0, 6.2); plt.ylim(0, 9.5)
plt.tight_layout()
plt.show()
Geometric corner points: [(np.float64(0.0), np.float64(0.0)), (np.float64(0.0), np.float64(6.0)), (np.float64(4.0), np.float64(0.0)
Objective values at corners: [np.float64(0.0), np.float64(30.0), np.float64(12.0), np.float64(27.0), np.float64(36.0)]
Best (geometric): x* = (np.float64(2.0), np.float64(6.0)) Z* = 36.0
SciPy (HiGHS) solution: x* = [2. 6.] Z* = 36.0
