[37] | 1 | function output = callmpt(interfacedata) |
---|
| 2 | |
---|
| 3 | % Author Johan Löfberg |
---|
| 4 | % $Id: callmpt.m,v 1.82 2006/09/12 14:08:33 joloef Exp $ |
---|
| 5 | |
---|
| 6 | % Speeds up solving LPs in mpmilp |
---|
| 7 | global mptOptions |
---|
| 8 | if ~isstruct(mptOptions) |
---|
| 9 | mpt_error |
---|
| 10 | end |
---|
| 11 | |
---|
| 12 | % Convert |
---|
| 13 | % interfacedata = pwa_linearize(interfacedata); |
---|
| 14 | Matrices = yalmip2mpt(interfacedata); |
---|
| 15 | |
---|
| 16 | % Get some MPT options |
---|
| 17 | options = interfacedata.options; |
---|
| 18 | options.mpt.lpsolver = mptOptions.lpsolver; |
---|
| 19 | options.mpt.milpsolver = mptOptions.milpsolver; |
---|
| 20 | options.mpt.verbose = options.verbose; |
---|
| 21 | |
---|
| 22 | if options.savedebug |
---|
| 23 | save mptdebug Matrices |
---|
| 24 | end |
---|
| 25 | |
---|
| 26 | if options.mp.unbounded |
---|
| 27 | Matrices = removeExplorationConstraints(Matrices); |
---|
| 28 | end |
---|
| 29 | |
---|
| 30 | if isempty(Matrices.binary_var_index) |
---|
| 31 | |
---|
| 32 | showprogress('Calling MPT',options.showprogress); |
---|
| 33 | solvertime = clock; |
---|
| 34 | if options.mp.presolve |
---|
| 35 | [Matrices.lb,Matrices.ub] = mpt_detect_and_improve_bounds(Matrices,Matrices.lb,Matrices.ub,Matrices.binary_var_index,options); |
---|
| 36 | end |
---|
| 37 | |
---|
| 38 | if any(Matrices.lb(end-Matrices.nx+1:end) == Matrices.ub(end-Matrices.nx+1:end)) |
---|
| 39 | model = []; |
---|
| 40 | else |
---|
| 41 | model = mpt_solvenode(Matrices,Matrices.lb,Matrices.ub,Matrices,[],options); |
---|
| 42 | end |
---|
| 43 | solvertime = etime(clock,solvertime); |
---|
| 44 | |
---|
| 45 | else |
---|
| 46 | % Pre-solve required on binary problems |
---|
| 47 | options.mp.presolve = 1; |
---|
| 48 | |
---|
| 49 | solvertime = clock; |
---|
| 50 | |
---|
| 51 | % if Matrices.qp & options.mp.algorithm == 3 |
---|
| 52 | % options.mp.algorithm = 1; |
---|
| 53 | % end |
---|
| 54 | % options.mp.algorithm = 3; |
---|
| 55 | |
---|
| 56 | switch options.mp.algorithm |
---|
| 57 | case 1 |
---|
| 58 | showprogress('Calling MPT via enumeration',options.showprogress); |
---|
| 59 | model = mpt_enumeration_mpmilp(Matrices,options); |
---|
| 60 | case 2 |
---|
| 61 | % Still experimental and just for fun. Not working! |
---|
| 62 | showprogress('Calling MPT via parametric B&B',options.showprogress); |
---|
| 63 | model = mpt_parbb(Matrices,options); |
---|
| 64 | |
---|
| 65 | case 3 |
---|
| 66 | showprogress('Calling MPT via delayed enumeration',options.showprogress); |
---|
| 67 | %Matrices = initialize_binary_equalities(Matrices) |
---|
| 68 | [Matrices.lb,Matrices.ub] = mpt_detect_and_improve_bounds(Matrices,Matrices.lb,Matrices.ub,Matrices.binary_var_index,options); |
---|
| 69 | model = mpt_de_mpmilp(Matrices,options,[]); |
---|
| 70 | |
---|
| 71 | otherwise |
---|
| 72 | end |
---|
| 73 | solvertime = etime(clock,solvertime); |
---|
| 74 | end |
---|
| 75 | |
---|
| 76 | if isempty(model) |
---|
| 77 | model = {model}; |
---|
| 78 | end |
---|
| 79 | |
---|
| 80 | if options.verbose |
---|
| 81 | if ~isempty(model{1}) |
---|
| 82 | if length(model) == 1 |
---|
| 83 | disp(['-> Generated 1 partition.']) |
---|
| 84 | else |
---|
| 85 | disp(['-> Generated ' num2str(length(model)) ' partitions.']) |
---|
| 86 | end |
---|
| 87 | end |
---|
| 88 | end |
---|
| 89 | |
---|
| 90 | problem = 0; |
---|
| 91 | infostr = yalmiperror(problem,'MPT'); |
---|
| 92 | |
---|
| 93 | % Save all data sent to solver? |
---|
| 94 | if options.savesolverinput |
---|
| 95 | solverinput.Matrices = Matrices; |
---|
| 96 | solverinput.options = []; |
---|
| 97 | else |
---|
| 98 | solverinput = []; |
---|
| 99 | end |
---|
| 100 | |
---|
| 101 | % Save all data from the solver? |
---|
| 102 | % This always done |
---|
| 103 | if options.savesolveroutput |
---|
| 104 | solveroutput.model = model; |
---|
| 105 | solveroutput.U = interfacedata.used_variables(Matrices.free_var);%(Matrices.free_var <= length( interfacedata.used_variables))); |
---|
| 106 | solveroutput.x = interfacedata.used_variables(Matrices.param_var); |
---|
| 107 | else |
---|
| 108 | solveroutput = []; |
---|
| 109 | end |
---|
| 110 | |
---|
| 111 | % Standard interface |
---|
| 112 | output.Primal = nan*ones(length(interfacedata.c),1); |
---|
| 113 | output.Dual = []; |
---|
| 114 | output.Slack = []; |
---|
| 115 | output.problem = problem; |
---|
| 116 | output.infostr = infostr; |
---|
| 117 | output.solverinput = solverinput; |
---|
| 118 | output.solveroutput= solveroutput; |
---|
| 119 | output.solvertime = solvertime; |
---|
| 120 | |
---|
| 121 | % |
---|
| 122 | % function M = pwa_linearize(M) |
---|
| 123 | % |
---|
| 124 | % if any(M.variabletype ~= 0) |
---|
| 125 | % [lb,ub] = findulb(M.F_struc,M.K); |
---|
| 126 | % M.lb = lb; |
---|
| 127 | % M.ub = ub; |
---|
| 128 | % nonlinear = find(M.variabletype ~= 0); |
---|
| 129 | % ok = 1; |
---|
| 130 | % for i = 1:length(nonlinear) |
---|
| 131 | % if nnz(M.monomtable(nonlinear(i),:)) == 1 |
---|
| 132 | % j = find(M.monomtable(nonlinear(i),:)); |
---|
| 133 | % if isinf(lb(j)) | isinf(ub(j)) |
---|
| 134 | % ok = 0; |
---|
| 135 | % end |
---|
| 136 | % else |
---|
| 137 | % ok = 0; |
---|
| 138 | % end |
---|
| 139 | % end |
---|
| 140 | % if ok |
---|
| 141 | % N = 5; |
---|
| 142 | % for i = 1:length(nonlinear) |
---|
| 143 | % j = find(M.monomtable(nonlinear(i),:)); |
---|
| 144 | % x = linspace(lb(j),ub(j),N); |
---|
| 145 | % y = x.^M.monomtable(nonlinear(i),j); |
---|
| 146 | % k = []; |
---|
| 147 | % m = []; |
---|
| 148 | % for r = 1:N-1 |
---|
| 149 | % km = polyfit(x(r:r+1),y(r:r+1),1); |
---|
| 150 | % k = [k km(1)]; |
---|
| 151 | % m = [m km(2)]; |
---|
| 152 | % end |
---|
| 153 | % % Redefine as linear |
---|
| 154 | % M.monomtable(nonlinear(i),:) = 0; |
---|
| 155 | % M.monomtable(nonlinear(i),nonlinear(i)) = 1; |
---|
| 156 | % M.variabletype(nonlinear(i)) = 0; |
---|
| 157 | % % Add new binaries |
---|
| 158 | % M.monomtable = blkdiag(M.monomtable,eye(N-1)); |
---|
| 159 | % M.variabletype = [M.variabletype zeros(1,N-1)]; |
---|
| 160 | % binaries = [length(M.monomtable)-N+2:length(M.monomtable)]; |
---|
| 161 | % M.binary_variables = [M.binary_variables binaries]; |
---|
| 162 | % % M.free_var = [M.free_var binaries]; |
---|
| 163 | % M.ub(nonlinear(i)) = max(y); |
---|
| 164 | % M.lb(nonlinear(i)) = min(y); |
---|
| 165 | % M.lb(binaries) = 0; |
---|
| 166 | % M.ub(binaries) = 1; |
---|
| 167 | % |
---|
| 168 | % M.F_struc = [zeros(1,size(M.F_struc,2));M.F_struc]; |
---|
| 169 | % M.F_struc(1,1) = 1; |
---|
| 170 | % M.F_struc(1,1+binaries) = -1; |
---|
| 171 | % M.K.f = M.K.f + 1; |
---|
| 172 | % A = []; |
---|
| 173 | % b = []; |
---|
| 174 | % bigM = 2*(max(y)-min(y)); |
---|
| 175 | % for r = 1:N-1 |
---|
| 176 | % A(end+1,binaries(r)) = -bigM; |
---|
| 177 | % A(end,nonlinear(i)) = -1; |
---|
| 178 | % A(end,j) = k(r); |
---|
| 179 | % b(end+1) = m(r)+bigM; |
---|
| 180 | % end |
---|
| 181 | % for r = 1:N-1 |
---|
| 182 | % A(end+1,binaries(r)) = -bigM; |
---|
| 183 | % A(end,nonlinear(i)) = 1; |
---|
| 184 | % A(end,j) = -k(r); |
---|
| 185 | % b(end+1) = bigM-m(r); |
---|
| 186 | % end |
---|
| 187 | % bigM = 20; |
---|
| 188 | % for r = 1:N-1 |
---|
| 189 | % A(end+1,binaries(r)) = -bigM; |
---|
| 190 | % A(end,j) = -1; |
---|
| 191 | % b(end+1) = bigM+x(r+1); |
---|
| 192 | % end |
---|
| 193 | % for r = 1:N-1 |
---|
| 194 | % A(end+1,binaries(r)) = -bigM; |
---|
| 195 | % A(end,j) = 1; |
---|
| 196 | % b(end+1) = bigM-x(r); |
---|
| 197 | % end |
---|
| 198 | % |
---|
| 199 | % M.F_struc = [M.F_struc(1:M.K.f,:);b(:) A;M.F_struc(M.K.f+1:end,:)]; |
---|
| 200 | % M.K.l = M.K.l + 1; |
---|
| 201 | % end |
---|
| 202 | % M.c(length(M.variabletype)) = 0; |
---|
| 203 | % M.Q(length(M.variabletype),length(M.variabletype)) = 0; |
---|
| 204 | % end |
---|
| 205 | % end |
---|
| 206 | % |
---|
| 207 | % |
---|
| 208 | % |
---|
| 209 | |
---|
| 210 | function Matrices = initialize_binary_equalities(Matrices) |
---|
| 211 | binary_var_index = Matrices.binary_var_index; |
---|
| 212 | notbinary_var_index = setdiff(1:Matrices.nu,binary_var_index); |
---|
| 213 | % Detect and extract pure binary equalities. Used for simple pruning |
---|
| 214 | nbin = length(binary_var_index); |
---|
| 215 | only_binary = ~any(Matrices.Aeq(:,notbinary_var_index),2); |
---|
| 216 | Matrices.Aeq_bin = Matrices.Aeq(find(only_binary),binary_var_index); |
---|
| 217 | Matrices.beq_bin = Matrices.beq(find(only_binary),:); |
---|
| 218 | |
---|
| 219 | function Matrices = removeExplorationConstraints(Matrices); |
---|
| 220 | candidates = find((~any(Matrices.G,2)) & (sum(Matrices.E | Matrices.E,2) == 1)); |
---|
| 221 | if ~isempty(candidates) |
---|
| 222 | Matrices.bndA = -Matrices.E(candidates,:); |
---|
| 223 | Matrices.bndb = Matrices.W(candidates,:); |
---|
| 224 | Matrices.G(candidates,:) = []; |
---|
| 225 | Matrices.E(candidates,:) = []; |
---|
| 226 | Matrices.W(candidates,:) = []; |
---|
| 227 | end |
---|