1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
|
function [x,z,info] = lin_prog(A,b,c);
% [x,z,info] = lin_prog(A,b,c);
%
% solves the linear program
%
% (primal) minimize c'*x
% subject to A*x + b >= 0
%
% (dual) maximize -b'*z
% subject to A'*z = c
% z >= 0
%
% Input arguments:
% A: nxm-matrix.
% b: n-vector.
% c: m-vector.
%
% Output arguments
% x: primal solution.
% z: dual solution.
% info: string. info = 'optimal', or 'may be infeasible', or
% 'may be unbounded below', or 'maxiters exceeded'.
%
%
% The code adds an upper bound to the primal constraints,
% and it really solves the problem
%
% (primal) minimize c'*x
% subject to A*x + b >= 0
% e'*(A*x + b) <= M
%
% (dual) maximize -b'*(z-w*e) - M*w
% subject to A'*(z-w*e) = c
% z >= 0, w >= 0
%
% with M = 1e3 * || [A b] ||_1
% and e is the vector with all components one.
%
% If no optimal solution can be found with strict inequality
% e'*(b-A*x) < M, then the original problem is infeasible, or
% unbounded below, or the default value of M is too small.
%
[n,m] = size(A);
if (size(b,1) < size(b,2)), b = b'; end;
if ( size(b,1) ~= n ) | (size(b,2) ~= 1),
error('b must be an n-vector.');
end;
if (size(c,1) < size(c,2)), c = c'; end;
if ( size(c,1) ~= m ) | (size(c,2) ~= 1),
error('c must be an m-vector.');
end;
M = 1e3 * max(sum(abs([A b])));
%
% phase 1
%
disp(' '); disp(' PHASE 1.');
abstol = 1e-8; nu = 10.0; maxiters = 100;
[x0,Z0,z0,ul,iters] = ...
phase1([b A], ones(n,1), M, nu, abstol, maxiters);
if iters == maxiters,
info = 'maxiters exceeded';
x = []; z = [];
return;
end;
if (ul(1) > 0.0)
info = 'may be infeasible';
x=[]; z=[];
return;
end;
%
% phase 2
%
disp(' '); disp(' PHASE 2.');
M = max(M, 1e3*sum(b+A*x0)); % M must be greater than e'(b-A*x0)
% for bigM.m
abstol = 1e-8; reltol = 1e-3; nu = 10.0; maxiters = 100;
[x,z,w,ul,iters] = ...
bigM([b A], ones(n,1), c, x0, M, nu, abstol, reltol, 0.0, maxiters);
if iters == maxiters
info = 'maxiters exceeded';
x=[]; z=[];
return;
end;
if sum(b+A*x) > 0.9*M
info = 'may be unbounded below';
x=[]; z=[];
return;
end;
info = 'optimal';
|
