This was my work in developing the plotregion function in Julia to show the feasible region of an LP.

I really do cheat though and just use the fact that the BFPs are the vertices.



In [4]:

    
# Pkg.add("QuickHull")
using QuickHull
using Combinatorics
using FixedSizeArrays









    



INFO: Recompiling stale cache file /home/juser/.julia/lib/v0.5/FixedSizeArrays.ji for module FixedSizeArrays.



In [10]:

    
collect(combinations(1:5,3))









    Out[10]:





10-element Array{Array{Int64,1},1}:
 [1,2,3]
 [1,2,4]
 [1,2,5]
 [1,3,4]
 [1,3,5]
 [1,4,5]
 [2,3,4]
 [2,3,5]
 [2,4,5]
 [3,4,5]



In [11]:

    
setdiff(1:5,[1,2,3])









    Out[11]:





2-element Array{Int64,1}:
 4
 5

%% Define the LP
% minimize    -x1 - 2*x2
% subject to  -2*x1 + x2 <= 2
%             -x1 + 2*x2 <= 7
%              x1        <= 3
%              x1, x2 >= 0
%
% Except, we need this in standard form:
% minimize f'*x subject to. A*x <= b

c = [-1 -2];
A = [-2 1; -1 2; 1 0];
b = [2; 7; 3];
lb = [0; 0];
ub = []; % No upper bound

%% Find the solution using linprog
x = linprog(c,A,b,[],[],lb);

%% Define the full LP with all the slacks.
c = [-1 -2 0 0 0];
A1 = [-2 1; -1 2; 1 0];
b = [2; 7; 3];
A = [A1 eye(3)];
%x = linprog(c,[],[],A,b,zeros(5,1));



In [26]:

    
function basic_feasible_point(A::Matrix,b::Vector,set::Vector)
    m,n = size(A)
    @assert length(set) == m "need more indices to define a BFP"
    binds = set # basic variable indices
    ninds = setdiff(1:size(A,1),binds) # non-basic
    B = A[:,binds]
    N = A[:,ninds]
    #cb = c[binds]
    #cn = c[ninds]
    
    if rank(B) != m
        return (:Infeasible, 0)
    end
    
    xb = B\b
    x = zeros(eltype(xb),n)
    x[binds] = xb
    x[ninds] = zero(eltype(xb))    
    
    if any(xb .< 0)
        return (:Infeasible, x)
    else
        #lam = B'\cb
        #sn = cn - N'*lam
        return (:Feasible, x)
    end
end
A1 = [-2.0 1; -1 2; 1 0]
b = [2.0; 7; 3]
A = [A1 eye(3)]
basic_feasible_point(A,b,[1,2,3])









    



WARNING: Method definition basic_feasible_point(Array{T<:Any, 2}, Array{T<:Any, 1}, Array{T<:Any, 1}) in module Main at In[25]:2 overwritten at In[26]:2.






    Out[26]:





(:Feasible,[3.0,5.0,3.0,0.0,0.0])



In [23]:

    
verts = Vector{Vec{2,Float64}}()
push!(verts, Vec(5.0,6.0))









    Out[23]:





1-element Array{FixedSizeArrays.Vec{2,Float64},1}:
 Vec(5.0,6.0)



In [33]:

    
using Plots









    



UndefVarError: ? not defined



In [40]:

    
"""
Plot the feasible polytope 
Ax = b, x >= 0
for the first two components of x.
"""
function plotregion(A::Matrix,b::Vector)
    m,n = size(A)
    verts = Vector{Vec{2,Float64}}()
    for inds in combinations(1:n,m)
        bfp=basic_feasible_point(A,b,inds)
        if bfp[1] == :Feasible
            push!(verts, Vec(bfp[2][1], bfp[2][2]) )
        end
    end
    hull = qhull(verts)
    plot(Shape(hull),fillalpha=0.5, fillcolor="grey", label="")
end
plotregion(A,b)









    



WARNING: Method definition plotregion(Array{T<:Any, 2}, Array{T<:Any, 1}) in module Main at In[39]:7 overwritten at In[40]:7.
WARNING: replacing docs for 'plotregion :: Tuple{Array{T,2},Array{T,1}}' in module 'Main'.






    Out[40]: