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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
| function [errmsg,Z,X,t,c,fail] = BNB20(fun,x0,xstat,xl,xu,a,b,aeq,beq,nonlc,setts,opts,varargin);
% BNB20 Finds the constrained minimum of a function of several possibly integer variables.
% Usage: [errmsg,Z,X,t,c,fail] =
% BNB20(fun,x0,xstatus,xlb,xub,A,B,Aeq,Beq,nonlcon,settings,options,P1,P2,...)
%
% BNB solves problems of the form:
% Minimize F(x) subject to: xlb <= x0 <=xub
% A*x <= B Aeq*x=Beq
% C(x)<=0 Ceq(x)=0
% x(i) is continuous for xstatus(i)=0
% x(i) integer for xstatus(i)= 1
% x(i) fixed for xstatus(i)=2
%
% BNB uses:
% Optimization Toolbox Version 2.0 (R11) 09-Oct-1998
% From this toolbox fmincon.m is called. For more info type help fmincon.
%
% fun is the function to be minimized and should return a scalar. F(x)=feval(fun,x).
% x0 is the starting point for x. x0 should be a column vector.
% xstatus is a column vector describing the status of every variable x(i).
% xlb and xub are column vectors with lower and upper bounds for x.
% A and Aeq are matrices for the linear constrains.
% B and Beq are column vectors for the linear constrains.
% nonlcon is the function for the nonlinear constrains.
% [C(x);Ceq(x)]=feval(nonlcon,x). Both C(x) and Ceq(x) should be column vectors.
%
% errmsg is a string containing an error message if BNB found an error in the input.
% Z is the scalar result of the minimization, X the values of the accompanying variables.
% t is the time elapsed while the algorithm BNB has run and c is the number of BNB cycles.
% fail is the number of nonconvergent leaf sub-problems.
%
% settings is a row vector with settings for BNB:
% settings(1) (standard 0) if 1: if the sub-problem does not converge do not branch it and
% raise fail by one. Normally BNB will always branch a nonconvergent sub-problem so it can
% try again to find a solution.
% A sub-problem that is a leaf of the branch-and-bound-tree cannot be branched. If such
% a problem does not converge it will be considered infeasible and fail will be raised by one.
% settings(2) is the handle of main BNB GUI. Leave empty.
% options is an options structure. For details type help optimset.
% options.maxSQPiter is a variable used by fmincon (if modified as described in bnb20.m).
% maxSQPiter cannot be set by optimset because it is not part of the standard options
% structure. maxSQPiter is 1000 by default.
% P1,P2,... are parameters to be passed to fun and nonlcon.
% F(x)=feval(fun,x,P1,P2,...). [C(x);Ceq(x)]=feval(nonlcon,x,P1,P2,...).
% Type edit BNB20 for more info.
% To get rid of bugs and to stop fmincon from hanging make the following chances:
%
% In optim/private/nlconst.m ($Revision: 1.20 $ $Date: 1998/08/24 13:46:15 $):
% Get EXITFLAG independent of verbosity.
% After the lines: disp(' less than 2*options.TolFun but constraints are not satisfied.')
% end
% EXITFLAG = -1;
% end
% end
% status=1;
% add the line: if (strncmp(howqp, 'i',1) & mg > 0), EXITFLAG = -1; end;
%
% In optim/private/qpsub.m ($Revision: 1.21 $ $Date: 1998/09/01 21:37:56 $):
% Stop qpsub from hanging.
% After the line: % Andy Grace 7-9-90. Mary Ann Branch 9-30-96.
% add the line: global maxSQPIter;
% and changed the line: maxSQPiters = Inf;
% to the line: if exist('maxSQPiter','var'), maxSQPiters = maxSQPiter; else maxSQPiters=inf; end;
% I guess there was a reason to put maxSQPiters at infinity, but this works fine for me.
global maxSQPiter;
% STEP 0 CHECKING INPUT
Z=[]; X=[]; t=0; c=0; fail=0;
if nargin<2, errmsg='BNB needs at least 2 input arguments.'; return; end;
if isempty(fun), errmsg='No fun found.'; return;
elseif ~ischar(fun), errmsg='fun must be a string.'; return; end;
if isempty(x0), errmsg='No x0 found.'; return;
elseif ~isnumeric(x0) | ~isreal(x0) | size(x0,2)>1
errmsg='x0 must be a real column vector.'; return;
end;
xstatus=zeros(size(x0));
if nargin>2 & ~isempty(xstat)
if isnumeric(xstat) & isreal(xstat) & all(size(xstat)<=size(x0))
if all(xstat==round(xstat) & 0<=xstat & xstat<=2)
xstatus(1:size(xstat))=xstat;
else errmsg='xstatus must consist of the integers 0,1 en 2.'; return; end;
else errmsg='xstatus must be a real column vector the same size as x0.'; return; end;
end;
xlb=zeros(size(x0));
xlb(find(xstatus==0))=-inf;
if nargin>3 & ~isempty(xl)
if isnumeric(xl) & isreal(xl) & all(size(xl)<=size(x0))
xlb(1:size(xl,1))=xl;
else errmsg='xlb must be a real column vector the same size as x0.'; return; end;
end;
if any(x0<xlb), errmsg='x0 must be in the range xlb <= x0.'; return;
elseif any(xstatus==1 & (~isfinite(xlb) | xlb~=round(xlb)))
errmsg='xlb(i) must be an integer if x(i) is an integer variabele.'; return;
end;
xlb(find(xstatus==2))=x0(find(xstatus==2));
xub=ones(size(x0));
xub(find(xstatus==0))=inf;
if nargin>4 & ~isempty(xu)
if isnumeric(xu) & isreal(xu) & all(size(xu)<=size(x0))
xub(1:size(xu,1))=xu;
else errmsg='xub must be a real column vector the same size as x0.'; return; end;
end;
if any(x0>xub), errmsg='x0 must be in the range x0 <=xub.'; return;
elseif any(xstatus==1 & (~isfinite(xub) | xub~=round(xub)))
errmsg='xub(i) must be an integer if x(i) is an integer variabale.'; return;
end;
xub(find(xstatus==2))=x0(find(xstatus==2));
A=[];
if nargin>5 & ~isempty(a)
if isnumeric(a) & isreal(a) & size(a,2)==size(x0,1), A=a;
else errmsg='Matrix A not correct.'; return; end;
end;
B=[];
if nargin>6 & ~isempty(b)
if isnumeric(b) & isreal(b) & all(size(b)==[size(A,1) 1]), B=b;
else errmsg='Column vector B not correct.'; return; end;
end;
if isempty(B) & ~isempty(A), B=zeros(size(A,1),1); end;
Aeq=[];
if nargin>7 & ~isempty(aeq)
if isnumeric(aeq) & isreal(aeq) & size(aeq,2)==size(x0,1), Aeq=aeq;
else errmsg='Matrix Aeq not correct.'; return; end;
end;
Beq=[];
if nargin>8 & ~isempty(beq)
if isnumeric(beq) & isreal(beq) & all(size(beq)==[size(Aeq,1) 1]), Beq=beq;
else errmsg='Column vector Beq not correct.'; return; end;
end;
if isempty(Beq) & ~isempty(Aeq), Beq=zeros(size(Aeq,1),1); end;
nonlcon='';
if nargin>9 & ~isempty(nonlc)
if ischar(nonlc), nonlcon=nonlc;
else errmsg='fun must be a string.'; return; end;
end;
settings = [0 0];
if nargin>10 & ~isempty(setts)
if isnumeric(setts) & isreal(setts) & all(size(setts)<=size(settings))
settings(setts~=0)=setts(setts~=0);
else errmsg='settings should be a row vector of length 1 or 2.'; return; end;
end;
maxSQPiter=1000;
options=optimset('fmincon');
if nargin>11 & ~isempty(opts)
if isstruct(opts)
if isfield(opts,'MaxSQPIter')
if isnumeric(opts.MaxSQPIter) & isreal(opts.MaxSQPIter) & ...
all(size(opts.MaxSQPIter)==1) & opts.MaxSQPIter>0 & ...
round(opts.MaxSQPIter)==opts.MaxSQPIter
maxSQPiter=opts.MaxSQPIter;
opts=rmfield(opts,'MaxSQPIter');
else errmsg='options.maxSQPiter must be an integer >0.'; return; end;
end;
options=optimset(options,opts);
else errmsg='options must be a structure.'; return; end;
end;
evalreturn=0;
eval(['z=',fun,'(x0,varargin{:});'],'errmsg=''fun caused error.''; evalreturn=1;');
if evalreturn==1, return; end;
if ~isempty(nonlcon)
eval(['[C, Ceq]=',nonlcon,'(x0,varargin{:});'],'errmsg=''nonlcon caused error.''; evalreturn=1;');
if evalreturn==1, return; end;
if size(C,2)>1 | size(Ceq,2)>1, errmsg='C en Ceq must be column vectors.'; return; end;
end;
% STEP 1 INITIALISATION
currentwarningstate=warning;
warning off;
tic;
lx = size(x0,1);
z_incumbent=inf;
x_incumbent=inf*ones(size(x0));
I = ceil(sum(log2(xub(find(xstatus==1))-xlb(find(xstatus==1))+1))+size(find(xstatus==1),1)+1);
stackx0=zeros(lx,I);
stackx0(:,1)=x0;
stackxlb=zeros(lx,I);
stackxlb(:,1)=xlb;
stackxub=zeros(lx,I);
stackxub(:,1)=xub;
stackdepth=zeros(1,I);
stackdepth(1,1)=1;
stacksize=1;
xchoice=zeros(size(x0));
if ~isempty(Aeq)
j=0;
for i=1:size(Aeq,1)
if Beq(i)==1 & all(Aeq(i,:)==0 | Aeq(i,:)==1)
J=find(Aeq(i,:)==1);
if all(xstatus(J)~=0 & xchoice(J)==0 & xlb(J)==0 & xub(J)==1)
if all(xstatus(J)~=2) | all(x0(J(find(xstatus(J)==2)))==0)
j=j+1;
xchoice(J)=j;
if sum(x0(J))==0, errmsg='x0 not correct.'; return; end;
end;
end;
end;
end;
end;
errx=optimget(options,'TolX');
handleupdate=[];
if ishandle(settings(2))
taghandlemain=get(settings(2),'Tag');
if strcmp(taghandlemain,'main BNB GUI')
handleupdate=guiupd;
handleupdatemsg=findobj(handleupdate,'Tag','updatemessage');
bnbguicb('hide main');
drawnow;
end;
end;
optionsdisplay=getfield(options,'Display');
if strcmp(optionsdisplay,'iter') | strcmp(optionsdisplay,'final')
show=1;
else show=0; end;
% STEP 2 TERMINIATION
while stacksize>0
c=c+1;
% STEP 3 LOADING OF CSP
x0=stackx0(:,stacksize);
xlb=stackxlb(:,stacksize);
xub=stackxub(:,stacksize);
x0(find(x0<xlb))=xlb(find(x0<xlb));
x0(find(x0>xub))=xub(find(x0>xub));
depth=stackdepth(1,stacksize);
stacksize=stacksize-1;
percdone=round(100*(1-sum(0.5.^(stackdepth(1:(stacksize+1))-1))));
% UPDATE FOR USER
if ishandle(handleupdate) & strcmp(get(handleupdate,'Tag'),'update BNB GUI')
t=toc;
updatemsg={ ...
sprintf('searched %3d %% of three',percdone) ...
sprintf('Z : %12.4e',z_incumbent) ...
sprintf('t : %12.1f secs',t) ...
sprintf('c : %12d cycles',c-1) ...
sprintf('fail : %12d cycles',fail)};
set(handleupdatemsg,'String',updatemsg);
drawnow;
else
disp(sprintf('*** searched %3d %% of three',percdone));
disp(sprintf('*** Z : %12.4e',z_incumbent));
disp(sprintf('*** t : %12.1f secs',t));
disp(sprintf('*** c : %12d cycles',c-1));
disp(sprintf('*** fail : %12d cycles',fail));
end;
% STEP 4 RELAXATION
[x z convflag]=fmincon(fun,x0,A,B,Aeq,Beq,xlb,xub,nonlcon,options,varargin{:});
% STEP 5 FATHOMING
K = find(xstatus==1 & xlb~=xub);
separation=1;
if convflag<0 | (convflag==0 & settings(1))
% FC 1
separation=0;
if show, disp('*** branch pruned'); end;
if convflag==0,
fail=fail+1;
if show, disp('*** not convergent'); end;
elseif show, disp('*** not feasible');
end;
elseif z>=z_incumbent & convflag>0
% FC 2
separation=0;
if show
disp('*** branch pruned');
disp('*** ghosted');
end;
elseif all(abs(round(x(K))-x(K))<errx) & convflag>0
% FC 3
z_incumbent = z;
x_incumbent = x;
separation = 0;
if show
disp('*** branch pruned');
disp('*** new best solution found');
end;
end;
% STEP 6 SELECTION
if separation == 1 & ~isempty(K)
dzsep=-1;
for i=1:size(K,1)
dxsepc = abs(round(x(K(i)))-x(K(i)));
if dxsepc>=errx | convflag==0
xsepc = x; xsepc(K(i))=round(x(K(i)));
dzsepc = abs(feval(fun,xsepc,varargin{:})-z);
if dzsepc>dzsep
dzsep=dzsepc;
ixsep=K(i);
end;
end;
end;
% STEP 7 SEPARATION
if xchoice(ixsep)==0
% XCHOICE==0
branch=1;
domain=[xlb(ixsep) xub(ixsep)];
sepdepth=depth;
while branch==1
xboundary=(domain(1)+domain(2))/2;
if x(ixsep)<xboundary
domainA=[domain(1) floor(xboundary)];
domainB=[floor(xboundary+1) domain(2)];
else
domainA=[floor(xboundary+1) domain(2)];
domainB=[domain(1) floor(xboundary)];
end;
sepdepth=sepdepth+1;
stacksize=stacksize+1;
stackx0(:,stacksize)=x;
stackxlb(:,stacksize)=xlb;
stackxlb(ixsep,stacksize)=domainB(1);
stackxub(:,stacksize)=xub;
stackxub(ixsep,stacksize)=domainB(2);
stackdepth(1,stacksize)=sepdepth;
if domainA(1)==domainA(2)
stacksize=stacksize+1;
stackx0(:,stacksize)=x;
stackxlb(:,stacksize)=xlb;
stackxlb(ixsep,stacksize)=domainA(1);
stackxub(:,stacksize)=xub;
stackxub(ixsep,stacksize)=domainA(2);
stackdepth(1,stacksize)=sepdepth;
branch=0;
else
domain=domainA;
branch=1;
end;
end;
else
% XCHOICE~=0
L=find(xchoice==xchoice(ixsep));
M=intersect(K,L);
[dummy,N]=sort(x(M));
part1=M(N(1:floor(size(N)/2))); part2=M(N(floor(size(N)/2)+1:size(N)));
sepdepth=depth+1;
stacksize=stacksize+1;
stackx0(:,stacksize)=x;
O = (1-sum(stackx0(part1,stacksize)))/size(part1,1);
stackx0(part1,stacksize)=stackx0(part1,stacksize)+O;
stackxlb(:,stacksize)=xlb;
stackxub(:,stacksize)=xub;
stackxub(part2,stacksize)=0;
stackdepth(1,stacksize)=sepdepth;
stacksize=stacksize+1;
stackx0(:,stacksize)=x;
O = (1-sum(stackx0(part2,stacksize)))/size(part2,1);
stackx0(part2,stacksize)=stackx0(part2,stacksize)+O;
stackxlb(:,stacksize)=xlb;
stackxub(:,stacksize)=xub;
stackxub(part1,stacksize)=0;
stackdepth(1,stacksize)=sepdepth;
end;
elseif separation==1 & isempty(K)
fail=fail+1;
if show
disp('*** branch pruned');
disp('*** leaf not convergent');
end;
end;
end;
% STEP 8 OUTPUT
t=toc;
Z = z_incumbent;
X = x_incumbent;
errmsg='';
if ishandle(handleupdate)
taghandleupdate=get(handleupdate,'Tag');
if strcmp(taghandleupdate,'update BNB GUI')
close(handleupdate);
end;
end;
eval(['warning ',currentwarningstate]);
********************************************************
normalement si je veux appliquer ce code à un model mathématique simple, je dois déclarer tt les parametres de fonction BNB20 (fun,x0,xstatus...ect).
avec fun c la fonction obective que je veux minimiser.
mais quand j essaie d ecrire :
fun= 4*x+3*y;
x0= [0 0];
xstatus= [1 1];
xlb=[0;0];
xub=[inf; inf];
A= [1;2;3;4];
B= [1;2]; |
executer une fonction
Répondre avec citation
Partager