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 |
---|