Linear Program

A little program to draw the polyhedron of constraints and the objective function:

For instance with the constraints:

$$ \frac{X_{1}}{4}+\frac{X_{2}}{3}\leq 90, \ \frac{X_{1}}{8}+\frac{X_{2}}{3}\leq 80, \ X_{1}\leq 200, \ X_{1},X_{2}\geq 0
$$
and the objective function to maximize: $7.8X_{1}+7.1X_{2}$

We type linprog([[1/4, 1/3, 90], [1/8, 1/3, 80], [1,0,200], [1,0,0], [0,1,0]], [7.8, 7.1], 300, 350, 'inf') and we get:

with:

Python
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Sep 15 12:28:06 2018
 
@author: moi
"""
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.animation as animation
 
fig, ax = plt.subplots()
 
def polyhed(fs,xmax, ymax, order):
    """
    Draws the polyhedron of constraints
    fs : the coefs of the constraints
    order : 'sup' to maximise, 'inf' to minimise
    """
    y = np.arange(0.0, ymax, 0.1)
    x = np.arange(0.0, xmax, 0.1)
    for f in fs:
        a, b, c = f
        if b == 0 and c != 0: # line x = c/a
            ax.fill_betweenx(y,c/a,xmax,  alpha=0.5)
        elif a*b != 0: # line y = -ax/b + c/b
            f = lambda u: -a*u/b + c/b
            yy = np.vectorize(f)(x)
            if order == "inf":
                ax.fill_between(x,yy,ymax,  alpha=0.5)
            if order == "sup":
                ax.fill_between(x,-ymax,yy,  alpha=0.5)
            ax.fill_between(x,yy,ymax,  where=yy >= 0, alpha=0.5)
            ax.plot(x,yy)
        elif a == 0 and b != 0: # line y = c/b
            ax.plot(x, [c/b]*10*xmax)
        elif a == 0: # line x = 0
            ax.plot(x,[0]*10*xmax)
 
def vert(fs, xmax, ymax):
    """
    Computes the vertices of the polyhedron
    """
    X, Y = [], []
    for i in range(len(fs)-1):
        for j in range(i+1, len(fs)):
            a, b, c = fs[i] # the two lines
            A, B, C = fs[j]
            if b!= 0 or B != 0: # both lines non trivial
                solx = (C*b - c*B)/(A*b - a*B)
                if b == 0: # first is a vertical line
                    soly = (-A*solx + C)/B
                else:  # first is non trivial
                    soly = (-a*solx + c)/b
                if solx <= xmax and soly <= ymax: # stay in the window
                    X.append(solx)
                    Y.append(soly)
                    ax.annotate('({:.1f},{:.1f})'.format(solx, soly), (solx, soly))
    ax.scatter(X, Y, s=80, facecolors='none', edgecolors='r') # draw the points
 
def linprog(fs, ob, xmax, ymax, order):
    """
    puts all together and drw an animation
    """
    polyhed(fs,xmax, ymax, order)
    vert(fs, xmax, ymax)
    x = np.arange(0.0, xmax, 0.01)
    ims = []
    axes = plt.gca()
    axes.set_xlim([0,xmax])
    axes.set_ylim([0,ymax])
    for pas in np.arange(64):
        u = pas*2*ymax*ob[1]/64
        def fob(t):
            a, b = ob
            return max(0, (-a*t + u)/b)
        ims.append(plt.plot(x, np.vectorize(fob)(x), linestyle='-.', linewidth=2))
    im_ani = animation.ArtistAnimation(fig, ims, interval=50, repeat_delay=3000, blit=True)
    im_ani.save('proglin.mp4', fps=4)
 
# Example :
linprog([[1/4, 1/3, 90], [1/8, 1/3, 80], [1,0,200], [1,0,0], [0,1,0]], [7.8, 7.1], 300, 350, 'inf')
linprog([[10,2,5], [1,2,1],[1,8,2],[1,0,0]],[13,8],2.75,1.1,"sup")
plt.show()

Tags

courtesy of webmatter.de