Optimization: principles and algorithms, by Michel Bierlaire
twoPhasesSimplex.m
Go to the documentation of this file.
1 %> \file
2 %> Algorithm 16.5: Two phases simplex algorithm with tableau to solve a linear optimization problem in standard from. Implementation of algorithm 16.5 of \cite Bier15-book
3 %>
4 %> @note Tested with \ref run1615.m
5 %> @note Tested with \ref run1616.m
6 %> @note Calls \ref simplexTableau
7 %> @note Calls \ref pivoting
8 %> @note Called by \ref transhipment
9 %> @note Called by \ref gomory
10 %>
11 %> @author Michel Bierlaire
12 %> @date Sun Mar 22 14:25:13 2015
13 %> @ingroup Algorithms
14 %> @ingroup chap16
15
16 %> Solve a linear optimization problem in standard form using the tableau simplex with two phases
17 %> \f[ \min_{x \in \mathbb{R}^n} c^T x
18 %> \f]
19 %> subject to
20 %> \f[ Ax = b \f]
21 %> and
22 %> \f[ x \geq 0 \f]
23 %> where \f$A\in\mathbb{R}^{m\times n}\f$, \f$b\in\mathbb{R}^m\f$ and
24 %> \f$c\in \mathbb{R}^n\f$.
25 %> @param A the constraint matrix
26 %> @param b the constraint right hand side
27 %> @param c the cost vector for the objective function
28 %> @return [x,copt,finaltableau]
29 %> @return x the optimal solution
30 %> @return copt the value of the objective function
31 %> @return finaltableau the optimal tableau
32 function [x,copt,finaltableau,feasible,bounded] = twoPhasesSimplex(A,b,c)
33
34  [m,n] = size(A) ;
35
36  if (m != size(b))
37  error("The number of rows in A and b do not match") ;
38  endif
39  if (n != size(c))
40  error("The number of columns in A and do not match the size of c") ;
41  endif
42
43 % Make sure that b >= 0
44  for i=1:m
45  if (b(i) < 0)
46  b(i) *= -1 ;
47  A(i,:) *= -1 ;
48  endif
49  endfor
50
51 % First tableau for Phase I
52  tab = [ A eye(m,m) b ; zeros(1,n+m+1) ];
53  for i = 1:n
54  tab(m+1,i) = -sum(tab(1:m,i));
55  end
56  tab(m+1,n+m+1) = -sum(tab(1:m,n+m+1)) ;
57
58 % The basic variables are variables n+1 to n+m
59  rowindex = [n+(1:m)] ;
60 % Solve the phase I problem
61  [opttableau,bounded,rowindex] = simplexTableau(tab,rowindex) ;
62
63  if (bounded == 0)
64  feasible = 1 ;
65  x=[] ;
66  copt = 0 ;
67  finaltableau = [] ;
68  bounded = 0 ;
69  return ;
70  endif
71  if (opttableau(m+1,n+m+1) > sqrt(eps))
72  printf("Optimal cost is %e > %e. No feasible solution exists.\n",opttableau(m+1,n+m+1),sqrt(eps)) ;
73  feasible = 0 ;
74  x=[] ;
75  copt = 0 ;
76  finaltableau = [] ;
77  bounded = 1 ;
78  return ;
79  endif
80  feasible = 1 ;
81  ok = 0 ;
82
83  rowstoberemoved = [] ;
84 % Remove the auxiliary variables from the basis
85  for i=n+1:n+m
86  basic = any(rowindex == i) ;
87  if (basic != 0)
88  rowpivot = find(opttableau(:,i)') ;
89  % Search the pivot
90  colpivot = -1 ;
91  k=1 ;
92  while (colpivot < 0 && k <= n)
93  if (abs(opttableau(rowpivot,k)) > eps)
94  colpivot = k ;
95  else
96  k = k + 1 ;
97  endif
98  end
99  if (colpivot < 0)
100  rowstoberemoved(end + 1) = rowpivot ;
101  else
102  opttableau = pivoting(opttableau,rowpivot,colpivot) ;
103  rowindex(rowpivot) = colpivot ;
104  endif
105  endif
106  end
107 % Remove the redundant constraints
108  for r = size(rowstoberemoved,2):-1:1
109  opttableau(rowstoberemoved(r),:) = [] ;
110  rowindex(rowstoberemoved(r)) = [] ;
111  endfor
112 % Remove the columns
113  finaltableau = opttableau ;
114  finaltableau(:,n+1:n+m) = [] ;
115
116  [m,n] = size(finaltableau) ;
117  m = m-1 ;
118  n = n-1 ;
119
120  cb = zeros(m,1) ;
121  for i=1:m
122  cb(i) = c(rowindex(i)) ;
123  end
124
125  for i=1:n
126  if (any(rowindex==i))
127  finaltableau(m+1,i) = 0 ;
128  else
129  finaltableau(m+1,i) = c(i) - cb' * finaltableau(1:m,i) ;
130  endif
131  end
132  finaltableau(m+1,n+1) = - cb' * finaltableau(1:m,n+1) ;
133  x = zeros(n,1) ;
134
135 % Phase II
136
137  [finaltableau,bounded,rowindex] = simplexTableau(finaltableau,rowindex) ;
138  x = zeros(n,1) ;
139  for j=1:m
140  x(rowindex(j)) = finaltableau(j,n+1) ;
141  endfor
142  copt = -finaltableau(m+1,n+1) ;
143 endfunction
144
function twoPhasesSimplex(in A, in b, in c)
Solve a linear optimization problem in standard form using the tableau simplex with two phases subje...
function simplexTableau(in tab, in rowindex)
Solve a linear optimization problem in standard form using the tableau simplex.
function gomory(in A, in b, in c)
Solve an integer optimization problem by generating valid inequalities using a Gomory cut: subject t...
function transhipment(in adj, in cost, in lb, in ub, in supply, in useGlpk)
Solve the transhipment problem with bound constraints.
function pivoting(in tab, in l, in j)
Pivot the tableau.
function simplex(in A, in b, in c, in basis)
Applies the simplex method to solve subject to and , where , , and .