#!/usr/bin/env python # -*- coding: utf-8 -*- # # Python script to illustrate the convergence of Karmarkar's algorithm on # a linear programming problem. # # http://en.wikipedia.org/wiki/Karmarkar%27s_algorithm # # Karmarkar's algorithm is an algorithm introduced by Narendra Karmarkar in 1984 # for solving linear programming problems. It was the first reasonably efficient # algorithm that solves these problems in polynomial time. # # Karmarkar's algorithm falls within the class of interior point methods: the # current guess for the solution does not follow the boundary of the feasible # set as in the simplex method, but it moves through the interior of the feasible # region, improving the approximation of the optimal solution by a definite # fraction with every iteration, and converging to an optimal solution with # rational data. # # Guillaume Jacquenot # 2015-05-03 # CC-BY-SA import numpy as np import matplotlib from matplotlib.pyplot import figure, show, rc, grid class ProblemInstance(): def __init__(self): n = 2 m = 11 self.A = np.zeros((m,n)) self.B = np.zeros((m,1)) self.C = np.array([[1],[1]]) self.A[:,1] = 1 for i in range(11): p = 0.1*i self.A[i,0] = 2.0*p self.B[i,0] = p*p + 1.0 class KarmarkarAlgorithm(): def __init__(self,A,B,C): self.maxIterations = 100 self.A = np.copy(A) self.B = np.copy(B) self.C = np.copy(C) self.n = len(C) self.m = len(B) self.AT = A.transpose() self.XT = None def isConvergeCriteronSatisfied(self, epsilon = 1e-8): if np.size(self.XT,1)<2: return False if np.linalg.norm(self.XT[:,-1]-self.XT[:,-2],2) < epsilon: return True def solve(self, X0=None): # No check is made for unbounded problem if X0 is None: X0 = np.zeros((self.n,1)) k = 0 X = np.copy(X0) self.XT = np.copy(X0) gamma = 0.5 for _ in range(self.maxIterations): if self.isConvergeCriteronSatisfied(): break V = self.B-np.dot(self.A,X) VM2 = np.linalg.matrix_power(np.diagflat(V),-2) hx = np.dot(np.linalg.matrix_power(np.dot(np.dot(self.AT,VM2),self.A),-1),self.C) hv = -np.dot(self.A,hx) coeff = np.infty for p in range(self.m): if hv[p,0]<0: coeff = np.min((coeff,-V[p,0]/hv[p,0])) alpha = gamma * coeff X += alpha*hx self.XT = np.concatenate((self.XT,X),axis=1) def makePlot(self,outputFilename = r'Karmarkar.svg', xs=-0.05, xe=+1.05): rc('grid', linewidth = 1, linestyle = '-', color = '#a0a0a0') rc('xtick', labelsize = 15) rc('ytick', labelsize = 15) rc('font',**{'family':'serif','serif':['Palatino'],'size':15}) rc('text', usetex=True) fig = figure() ax = fig.add_axes([0.12, 0.12, 0.76, 0.76]) grid(True) ylimMin = -0.05 ylimMax = +1.05 xsOri = xs xeOri = xe for i in range(np.size(self.A,0)): xs = xsOri xe = xeOri a = -self.A[i,0]/self.A[i,1] b = +self.B[i,0]/self.A[i,1] ys = a*xs+b ye = a*xe+b if ys>ylimMax: ys = ylimMax xs = (ylimMax-b)/a if ye<ylimMin: ye = ylimMin xe = (ylimMin-b)/a ax.plot([xs,xe], [ys,ye], \ lw = 1, ls = '--', color = 'b') ax.set_xlim((xs,xe)) ax.plot(self.XT[0,:], self.XT[1,:], \ lw = 1, ls = '-', color = 'r', marker = '.') ax.plot(self.XT[0,-1], self.XT[1,-1], \ lw = 1, ls = '-', color = 'r', marker = 'o') ax.set_xlim((ylimMin,ylimMax)) ax.set_ylim((ylimMin,ylimMax)) ax.set_aspect('equal') ax.set_xlabel('$x_1$',rotation = 0) ax.set_ylabel('$x_2$',rotation = 0) ax.set_title(r'$\max x_1+x_2\textrm{ s.t. }2px_1+x_2\le p^2+1\textrm{, }\forall p \in [0.0,0.1,...,1.0]$', fontsize=15) fig.savefig(outputFilename) fig.show() if __name__ == "__main__": p = ProblemInstance() k = KarmarkarAlgorithm(p.A,p.B,p.C) k.solve(X0 = np.zeros((2,1))) k.makePlot(outputFilename = r'Karmarkar.svg', xs=-0.05, xe=+1.05)