標準形の2次計画問題
minimize z = 0.5x'Qx + c'x
subject to Ax = b, x >= 0
を解くためのOctave関数を紹介します．
octave:1> A = [1 2 1 0; > 1 1 0 1] A = 1 2 1 0 1 1 0 1 octave:2> b = [7;5] b = 7 5 octave:3> Q = [1 0.5 0 0; > 0.5 1 0 0; > 0 0 0 0; > 0 0 0 0] Q = 1.00000 0.50000 0.00000 0.00000 0.50000 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 octave:4> c = -[10;11;0;0] c = -10 -11 -0 -0 octave:5> [x,zmin,is_bounded,sol] = qp(Q,c,A,b) x = 3 2 0 0 zmin = -42.500 is_bounded = 1 sol = 1
function [optx,zmin,is_bounded,sol] = qp(Q,c,A,b)
% [optx,zmin,is_bounded,sol] = qp(Q,c,A,b)
% This program solves the standard quadratic programming problem:
% minimize z = 0.5x'Qx + c'x
% subject to Ax = b, x >= 0
% using the algorithm of Dantzig and Van de Panne, a simplex-like method.
% optx: an optimal solution.
% zmin: the optimal value.
% is_bounded: 1 if the solution is bounded; 0 if unbounded.
% sol: 1 if the problem is solvable; 0 if unsolvable.
%
% Actually the program solves the equivalent problem:
% maximize z = -0.5x'Cx + p'x
% subject to Ax = b, x >= 0
% where C = Q, p = -c.
C = Q;
p = -c;
[n,m] = size(A);
c1 = zeros(m,1);
[optx,zmin,is_bounded,sol,basis] = lp(c1,A,b); % getting a feasible solution
if sol != 0
basis = sort(basis); % indices of basis variables
% indices of nonbasis variables
nonbasis = 1:m;
for i=basis
nonbasis(i) = 0;
endfor
for i=m:-1:1
if nonbasis(i) == 0
nonbasis(i) = [];
endif
endfor
index = [basis nonbasis];
% rearrangement of the matrices according to the indices
for j = 1:m
A0(:,j) = A(:,index(j));
p0(j) = p(index(j));
Ctemp(:,j) = C(:,index(j));
endfor
C0 = zeros(m,m);
for i = 1:m
C0(i,:) = Ctemp(index(i),:);
endfor
% making the initial simplex-like table
H = [eye(m,m) -A0' -C0 zeros(m,n) -p0';
zeros(n,m) zeros(n,n) A0 eye(n,n) b];
for ip = 1:n
jp = m+ip;
H = pivot(H,ip,jp);
endfor
for i = 1:n
ip = m+i;
jp = m+n+i;
H = pivot(H,ip,jp);
endfor
H(1:n,:) = [];
H(:,m+1:m+n) = [];
H(:,m+m+1:m+m+n) = [];
n1 = m-n;
m2 = m+m;
m3 = m2+1;
basis = n+1:m;
basis = [basis [m+1:m+n]];
us = zeros(1,m);
xs = us;
for i=basis
if i > m
xs(i-m) = 1;
else
us(i) = 1;
endif
endfor
opt = 0;
while opt == 0
if us*xs' == 0 % in standard case
umin = 1;
ns = 0;
for i = 1:m
if basis(i) <= m
if H(i,m3) < umin
umin = H(i,m3);
ns = i;
endif
endif
endfor
if umin >= 0
%disp('opt')
opt = 1; % an optimal solution is obtained
else
k = basis(ns)+m; % k:index to be in the basis. basis(ns):index of candidate to be out the basis
endif
else % in nonstandard case
uxsum = us + xs;
for i = 1:m
if uxsum(i) == 0
k = i;
break;
endif
endfor
for i = 1:m
if basis(i) <= m
if uxsum(basis(i)) == 2
ns = i;
break;
endif
endif
endfor
end
if opt == 0
% finding index to be out the basis
amin = Inf;
i0 = 0;
if H(ns,k) != 0
a = H(ns,m3)/H(ns,k);
if (a >= 0) & (a < amin)
amin = a;
i0 = ns;
endif
endif
for i = 1:m
if basis(i) > m
if H(i,k) != 0
a = H(i,m3)/H(i,k);
if (a >=0) & (a < amin)
amin = a;
i0 = i;
endif
endif
endif
endfor
% letting basis(i0) be out the basis and k in
if i0 == 0
%disp('unbounded')
opt = 1;
is_bounded = 0; % unbounded solution
else
basis(i0) = k;
H = pivot(H,i0,k);
us = zeros(1,m);
xs = us;
for i=basis
if i > m
xs(i-m) = 1;
else
us(i) = 1;
endif
endfor
endif
endif
endwhile
u = zeros(m,1);
x = u;
for i = 1:m
if basis(i) > m
x(basis(i)-m) = H(i,m3);
else
u(basis(i)) = H(i,m3);
endif
endfor
for i = 1:m
optx(index(i)) = x(i);
optu(index(i)) = u(i);
endfor
zmin = 0.5*optx'*Q*optx + c'*optx; % optimal value of z
endif
endfunction
octave:1> A = [1 2 1 0;
> 1 1 0 1]
A =
1 2 1 0
1 1 0 1
octave:2> b = [7;5]
b =
7
5
octave:3> Q = [1 0.5 0 0;
> 0.5 1 0 0;
> 0 0 0 0;
> 0 0 0 0]
Q =
1.00000 0.50000 0.00000 0.00000
0.50000 1.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000
octave:4> c = [-10;-11;0;0]
c =
-10
-11
0
0
octave:5> [x,y,obhis] = solqp(Q,A,b,c)
Search for a feasible point:
gap = 0
Search for an optimal solution:
obvalue = -31.988
obvalue = -41.816
lower = -42.863
obvalue = -42.468
obvalue = -42.495
lower = -42.505
obvalue = -42.500
obvalue = -42.500
A (local) optimal solution is found.
x =
3.0000e+00
2.0000e+00
1.5917e-06
8.5051e-06
y =
-1.5001
-4.4998
obhis =
-31.988 -41.816 -42.468 -42.495 -42.500 -42.500
octave:6> zmin = 0.5*x'*Q*x + c'*x
zmin = -42.500