In [19]:

    

using Convex
using SCS
using Plots


    

In [2]:

    

include("plotregion.jl")


    

    
Out[2]:





PlotRegion

In [3]:

    

A1 = [-2.0 1; 
      -1 2; 
       1 0]
b = [2.0; 7; 3]

A = [A1 eye(3)] # form the problem with slacks.

PlotRegion.plotregion(A,b)


    

    





    







    
Out[3]:




    

    

In [6]:

    

# Convert the problem into standard form.
cs = [-1 -2 0 0 0]'
AS = A


    

    
Out[6]:





3×5 Array{Float64,2}:
 -2.0  1.0  1.0  0.0  0.0
 -1.0  2.0  0.0  1.0  0.0
  1.0  0.0  0.0  0.0  1.0

In [37]:

    

""" Solve the central-path problem for interior point methods. """
function ip_central(c,A,b,tau)
    x = Variable(length(cs))
    p = minimize(c'*x - tau*sum(log(x)))
    p.constraints += A*x == b
    p.constraints += x >= 0
    solve!(p, SCSSolver(verbose=false))
    return x.value, p
end
ip_central(cs,AS,b,10.0)[1]


    

    




WARNING: Method definition ip_central(Any, Any, Any, Any) in module Main at In[13]:3 overwritten at In[37]:3.
WARNING: replacing docs for 'ip_central :: Tuple{Any,Any,Any,Any}' in module 'Main'.






    
Out[37]:





5×1 Array{Float64,2}:
 2.19431 
 2.19697 
 4.19166 
 4.80038 
 0.805689

In [26]:

    

A1 = [-2.0 1; 
      -1 2; 
       1 0]
b = [2.0; 7; 3]

AS = [A1 eye(3)] # form the problem with slacks.
cs = [-1 -2 0 0 0]'

taus = vec([10 7.5 5 3.5 2 1 0.75 0.5 0.35 0.20 logspace(-1,-7,10)'])

p = PlotRegion.plotregion(AS,b)
for tau in taus
    x = ip_central(cs, AS, b, tau)[1]
    scatter!([x[1]],[x[2]],label="", color="red")
end
p


    

    




WARNING: Problem status UnknownError; solution may be inaccurate.
WARNING: Problem status UnknownError; solution may be inaccurate.
WARNING: Problem status UnknownError; solution may be inaccurate.
WARNING: Problem status UnknownError; solution may be inaccurate.
WARNING: Problem status UnknownError; solution may be inaccurate.
WARNING: Problem status UnknownError; solution may be inaccurate.






    
Out[26]:




    

    

In [67]:

    

tau = 1.0
x0,prob = ip_central(cs,AS,b,tau)
x = copy(x0)
x[1] += 0.1
x[2] += 0.1

lam = vec(prob.constraints[1].dual)
@show lam

# show the region and the starting point
plt = PlotRegion.plotregion(AS,b)
scatter!([x[1]],[x[2]],label="", color="red")

# compute the steps 
n = length(cs)
m = size(AS,1)

s = tau./x
J = [zeros(n,n) AS'         eye(n); 
     AS        zeros(m,m)   zeros(m,n);
    diagm(vec(s))  zeros(n,m) diagm(vec(x))]

mu = dot(x,s)/n
sigma = 0.5

F = [s + AS'*lam - cs; AS*x - b; x.*s]

Fc = [s + AS'*lam - cs; AS*x - b; x.*s - sigma*mu ]

p = J\-F
pc = J\-Fc

plot!([x[1];x[1] + p[1]], [x[2];x[2] + p[2]], label="Affine")
plot!([x[1];x[1] + pc[1]], [x[2];x[2] + pc[2]], label = "Centered")

plt


    

    




lam = [-0.330829,-0.950675,-2.98753]






    
Out[67]:




    

    

In [65]:

    

xf = [x; lam; s];

[xf+p xf+pc]


    

    
Out[65]:





13×2 Array{Float64,2}:
  2.88968     2.77744 
  4.95807     4.63237 
  2.82129     2.92251 
 -0.0264644   0.512707
  0.110321    0.222562
 -0.0221681  -0.176443
 -0.97472    -0.962701
 -2.00279    -2.49481 
 -0.0162647   0.179225
 -0.0283913   0.101844
  0.0221681   0.176443
  0.97472     0.962701
  2.00279     2.49481 

In [63]:

    

taus = vec([10 7.5 5 3.5 2 1 0.75 0.5 0.35 0.20 logspace(-1,-7,10)'])

#p = PlotRegion.plotregion(AS,b)
for tau in taus
    x = ip_central(cs, AS, b, tau)[1]
    scatter!([x[1]],[x[2]],label="", color="red")
end
plt


    

    




WARNING: Problem status UnknownError; solution may be inaccurate.
WARNING: Problem status UnknownError; solution may be inaccurate.
WARNING: Problem status UnknownError; solution may be inaccurate.
WARNING: Problem status UnknownError; solution may be inaccurate.
WARNING: Problem status UnknownError; solution may be inaccurate.
WARNING: Problem status UnknownError; solution may be inaccurate.






    
Out[63]: