Dawn@4
|
1 function [x,fval,exitflag,output]=fminsearchbnd(fun,x0,LB,UB,options,varargin)
|
Dawn@4
|
2 % FMINSEARCHBND: FMINSEARCH, but with bound constraints by transformation
|
Dawn@4
|
3 % usage: x=FMINSEARCHBND(fun,x0)
|
Dawn@4
|
4 % usage: x=FMINSEARCHBND(fun,x0,LB)
|
Dawn@4
|
5 % usage: x=FMINSEARCHBND(fun,x0,LB,UB)
|
Dawn@4
|
6 % usage: x=FMINSEARCHBND(fun,x0,LB,UB,options)
|
Dawn@4
|
7 % usage: x=FMINSEARCHBND(fun,x0,LB,UB,options,p1,p2,...)
|
Dawn@4
|
8 % usage: [x,fval,exitflag,output]=FMINSEARCHBND(fun,x0,...)
|
Dawn@4
|
9 %
|
Dawn@4
|
10 % arguments:
|
Dawn@4
|
11 % fun, x0, options - see the help for FMINSEARCH
|
Dawn@4
|
12 %
|
Dawn@4
|
13 % LB - lower bound vector or array, must be the same size as x0
|
Dawn@4
|
14 %
|
Dawn@4
|
15 % If no lower bounds exist for one of the variables, then
|
Dawn@4
|
16 % supply -inf for that variable.
|
Dawn@4
|
17 %
|
Dawn@4
|
18 % If no lower bounds at all, then LB may be left empty.
|
Dawn@4
|
19 %
|
Dawn@4
|
20 % Variables may be fixed in value by setting the corresponding
|
Dawn@4
|
21 % lower and upper bounds to exactly the same value.
|
Dawn@4
|
22 %
|
Dawn@4
|
23 % UB - upper bound vector or array, must be the same size as x0
|
Dawn@4
|
24 %
|
Dawn@4
|
25 % If no upper bounds exist for one of the variables, then
|
Dawn@4
|
26 % supply +inf for that variable.
|
Dawn@4
|
27 %
|
Dawn@4
|
28 % If no upper bounds at all, then UB may be left empty.
|
Dawn@4
|
29 %
|
Dawn@4
|
30 % Variables may be fixed in value by setting the corresponding
|
Dawn@4
|
31 % lower and upper bounds to exactly the same value.
|
Dawn@4
|
32 %
|
Dawn@4
|
33 % Notes:
|
Dawn@4
|
34 %
|
Dawn@4
|
35 % If options is supplied, then TolX will apply to the transformed
|
Dawn@4
|
36 % variables. All other FMINSEARCH parameters should be unaffected.
|
Dawn@4
|
37 %
|
Dawn@4
|
38 % Variables which are constrained by both a lower and an upper
|
Dawn@4
|
39 % bound will use a sin transformation. Those constrained by
|
Dawn@4
|
40 % only a lower or an upper bound will use a quadratic
|
Dawn@4
|
41 % transformation, and unconstrained variables will be left alone.
|
Dawn@4
|
42 %
|
Dawn@4
|
43 % Variables may be fixed by setting their respective bounds equal.
|
Dawn@4
|
44 % In this case, the problem will be reduced in size for FMINSEARCH.
|
Dawn@4
|
45 %
|
Dawn@4
|
46 % The bounds are inclusive inequalities, which admit the
|
Dawn@4
|
47 % boundary values themselves, but will not permit ANY function
|
Dawn@4
|
48 % evaluations outside the bounds. These constraints are strictly
|
Dawn@4
|
49 % followed.
|
Dawn@4
|
50 %
|
Dawn@4
|
51 % If your problem has an EXCLUSIVE (strict) constraint which will
|
Dawn@4
|
52 % not admit evaluation at the bound itself, then you must provide
|
Dawn@4
|
53 % a slightly offset bound. An example of this is a function which
|
Dawn@4
|
54 % contains the log of one of its parameters. If you constrain the
|
Dawn@4
|
55 % variable to have a lower bound of zero, then FMINSEARCHBND may
|
Dawn@4
|
56 % try to evaluate the function exactly at zero.
|
Dawn@4
|
57 %
|
Dawn@4
|
58 %
|
Dawn@4
|
59 % Example usage:
|
Dawn@4
|
60 % rosen = @(x) (1-x(1)).^2 + 105*(x(2)-x(1).^2).^2;
|
Dawn@4
|
61 %
|
Dawn@4
|
62 % fminsearch(rosen,[3 3]) % unconstrained
|
Dawn@4
|
63 % ans =
|
Dawn@4
|
64 % 1.0000 1.0000
|
Dawn@4
|
65 %
|
Dawn@4
|
66 % fminsearchbnd(rosen,[3 3],[2 2],[]) % constrained
|
Dawn@4
|
67 % ans =
|
Dawn@4
|
68 % 2.0000 4.0000
|
Dawn@4
|
69 %
|
Dawn@4
|
70 % See test_main.m for other examples of use.
|
Dawn@4
|
71 %
|
Dawn@4
|
72 %
|
Dawn@4
|
73 % See also: fminsearch, fminspleas
|
Dawn@4
|
74 %
|
Dawn@4
|
75 %
|
Dawn@4
|
76 % Author: John D'Errico
|
Dawn@4
|
77 % E-mail: woodchips@rochester.rr.com
|
Dawn@4
|
78 % Release: 4
|
Dawn@4
|
79 % Release date: 7/23/06
|
Dawn@4
|
80
|
Dawn@4
|
81 % size checks
|
Dawn@4
|
82 xsize = size(x0);
|
Dawn@4
|
83 x0 = x0(:);
|
Dawn@4
|
84 n=length(x0);
|
Dawn@4
|
85
|
Dawn@4
|
86 if (nargin<3) || isempty(LB)
|
Dawn@4
|
87 LB = repmat(-inf,n,1);
|
Dawn@4
|
88 else
|
Dawn@4
|
89 LB = LB(:);
|
Dawn@4
|
90 end
|
Dawn@4
|
91 if (nargin<4) || isempty(UB)
|
Dawn@4
|
92 UB = repmat(inf,n,1);
|
Dawn@4
|
93 else
|
Dawn@4
|
94 UB = UB(:);
|
Dawn@4
|
95 end
|
Dawn@4
|
96
|
Dawn@4
|
97 if (n~=length(LB)) || (n~=length(UB))
|
Dawn@4
|
98 error 'x0 is incompatible in size with either LB or UB.'
|
Dawn@4
|
99 end
|
Dawn@4
|
100
|
Dawn@4
|
101 % set default options if necessary
|
Dawn@4
|
102 if (nargin<5) || isempty(options)
|
Dawn@4
|
103 options = optimset('fminsearch');
|
Dawn@4
|
104 end
|
Dawn@4
|
105
|
Dawn@4
|
106 % stuff into a struct to pass around
|
Dawn@4
|
107 params.args = varargin;
|
Dawn@4
|
108 params.LB = LB;
|
Dawn@4
|
109 params.UB = UB;
|
Dawn@4
|
110 params.fun = fun;
|
Dawn@4
|
111 params.n = n;
|
Dawn@4
|
112 params.OutputFcn = [];
|
Dawn@4
|
113
|
Dawn@4
|
114 % 0 --> unconstrained variable
|
Dawn@4
|
115 % 1 --> lower bound only
|
Dawn@4
|
116 % 2 --> upper bound only
|
Dawn@4
|
117 % 3 --> dual finite bounds
|
Dawn@4
|
118 % 4 --> fixed variable
|
Dawn@4
|
119 params.BoundClass = zeros(n,1);
|
Dawn@4
|
120 for i=1:n
|
Dawn@4
|
121 k = isfinite(LB(i)) + 2*isfinite(UB(i));
|
Dawn@4
|
122 params.BoundClass(i) = k;
|
Dawn@4
|
123 if (k==3) && (LB(i)==UB(i))
|
Dawn@4
|
124 params.BoundClass(i) = 4;
|
Dawn@4
|
125 end
|
Dawn@4
|
126 end
|
Dawn@4
|
127
|
Dawn@4
|
128 % transform starting values into their unconstrained
|
Dawn@4
|
129 % surrogates. Check for infeasible starting guesses.
|
Dawn@4
|
130 x0u = x0;
|
Dawn@4
|
131 k=1;
|
Dawn@4
|
132 for i = 1:n
|
Dawn@4
|
133 switch params.BoundClass(i)
|
Dawn@4
|
134 case 1
|
Dawn@4
|
135 % lower bound only
|
Dawn@4
|
136 if x0(i)<=LB(i)
|
Dawn@4
|
137 % infeasible starting value. Use bound.
|
Dawn@4
|
138 x0u(k) = 0;
|
Dawn@4
|
139 else
|
Dawn@4
|
140 x0u(k) = sqrt(x0(i) - LB(i));
|
Dawn@4
|
141 end
|
Dawn@4
|
142
|
Dawn@4
|
143 % increment k
|
Dawn@4
|
144 k=k+1;
|
Dawn@4
|
145 case 2
|
Dawn@4
|
146 % upper bound only
|
Dawn@4
|
147 if x0(i)>=UB(i)
|
Dawn@4
|
148 % infeasible starting value. use bound.
|
Dawn@4
|
149 x0u(k) = 0;
|
Dawn@4
|
150 else
|
Dawn@4
|
151 x0u(k) = sqrt(UB(i) - x0(i));
|
Dawn@4
|
152 end
|
Dawn@4
|
153
|
Dawn@4
|
154 % increment k
|
Dawn@4
|
155 k=k+1;
|
Dawn@4
|
156 case 3
|
Dawn@4
|
157 % lower and upper bounds
|
Dawn@4
|
158 if x0(i)<=LB(i)
|
Dawn@4
|
159 % infeasible starting value
|
Dawn@4
|
160 x0u(k) = -pi/2;
|
Dawn@4
|
161 elseif x0(i)>=UB(i)
|
Dawn@4
|
162 % infeasible starting value
|
Dawn@4
|
163 x0u(k) = pi/2;
|
Dawn@4
|
164 else
|
Dawn@4
|
165 x0u(k) = 2*(x0(i) - LB(i))/(UB(i)-LB(i)) - 1;
|
Dawn@4
|
166 % shift by 2*pi to avoid problems at zero in fminsearch
|
Dawn@4
|
167 % otherwise, the initial simplex is vanishingly small
|
Dawn@4
|
168 x0u(k) = 2*pi+asin(max(-1,min(1,x0u(k))));
|
Dawn@4
|
169 end
|
Dawn@4
|
170
|
Dawn@4
|
171 % increment k
|
Dawn@4
|
172 k=k+1;
|
Dawn@4
|
173 case 0
|
Dawn@4
|
174 % unconstrained variable. x0u(i) is set.
|
Dawn@4
|
175 x0u(k) = x0(i);
|
Dawn@4
|
176
|
Dawn@4
|
177 % increment k
|
Dawn@4
|
178 k=k+1;
|
Dawn@4
|
179 case 4
|
Dawn@4
|
180 % fixed variable. drop it before fminsearch sees it.
|
Dawn@4
|
181 % k is not incremented for this variable.
|
Dawn@4
|
182 end
|
Dawn@4
|
183
|
Dawn@4
|
184 end
|
Dawn@4
|
185 % if any of the unknowns were fixed, then we need to shorten
|
Dawn@4
|
186 % x0u now.
|
Dawn@4
|
187 if k<=n
|
Dawn@4
|
188 x0u(k:n) = [];
|
Dawn@4
|
189 end
|
Dawn@4
|
190
|
Dawn@4
|
191 % were all the variables fixed?
|
Dawn@4
|
192 if isempty(x0u)
|
Dawn@4
|
193 % All variables were fixed. quit immediately, setting the
|
Dawn@4
|
194 % appropriate parameters, then return.
|
Dawn@4
|
195
|
Dawn@4
|
196 % undo the variable transformations into the original space
|
Dawn@4
|
197 x = xtransform(x0u,params);
|
Dawn@4
|
198
|
Dawn@4
|
199 % final reshape
|
Dawn@4
|
200 x = reshape(x,xsize);
|
Dawn@4
|
201
|
Dawn@4
|
202 % stuff fval with the final value
|
Dawn@4
|
203 fval = feval(params.fun,x,params.args{:});
|
Dawn@4
|
204
|
Dawn@4
|
205 % fminsearchbnd was not called
|
Dawn@4
|
206 exitflag = 0;
|
Dawn@4
|
207
|
Dawn@4
|
208 output.iterations = 0;
|
Dawn@4
|
209 output.funcount = 1;
|
Dawn@4
|
210 output.algorithm = 'fminsearch';
|
Dawn@4
|
211 output.message = 'All variables were held fixed by the applied bounds';
|
Dawn@4
|
212
|
Dawn@4
|
213 % return with no call at all to fminsearch
|
Dawn@4
|
214 return
|
Dawn@4
|
215 end
|
Dawn@4
|
216
|
Dawn@4
|
217 % Check for an outputfcn. If there is any, then substitute my
|
Dawn@4
|
218 % own wrapper function.
|
Dawn@4
|
219 if ~isempty(options.OutputFcn)
|
Dawn@4
|
220 params.OutputFcn = options.OutputFcn;
|
Dawn@4
|
221 options.OutputFcn = @outfun_wrapper;
|
Dawn@4
|
222 end
|
Dawn@4
|
223
|
Dawn@4
|
224 % now we can call fminsearch, but with our own
|
Dawn@4
|
225 % intra-objective function.
|
Dawn@4
|
226 [xu,fval,exitflag,output] = fminsearch(@intrafun,x0u,options,params);
|
Dawn@4
|
227
|
Dawn@4
|
228 % undo the variable transformations into the original space
|
Dawn@4
|
229 x = xtransform(xu,params);
|
Dawn@4
|
230
|
Dawn@4
|
231 % final reshape
|
Dawn@4
|
232 x = reshape(x,xsize);
|
Dawn@4
|
233
|
Dawn@4
|
234 % Use a nested function as the OutputFcn wrapper
|
Dawn@4
|
235 function stop = outfun_wrapper(x,varargin);
|
Dawn@4
|
236 % we need to transform x first
|
Dawn@4
|
237 xtrans = xtransform(x,params);
|
Dawn@4
|
238
|
Dawn@4
|
239 % then call the user supplied OutputFcn
|
Dawn@4
|
240 stop = params.OutputFcn(xtrans,varargin{1:(end-1)});
|
Dawn@4
|
241
|
Dawn@4
|
242 end
|
Dawn@4
|
243
|
Dawn@4
|
244 end % mainline end
|
Dawn@4
|
245
|
Dawn@4
|
246 % ======================================
|
Dawn@4
|
247 % ========= begin subfunctions =========
|
Dawn@4
|
248 % ======================================
|
Dawn@4
|
249 function fval = intrafun(x,params)
|
Dawn@4
|
250 % transform variables, then call original function
|
Dawn@4
|
251
|
Dawn@4
|
252 % transform
|
Dawn@4
|
253 xtrans = xtransform(x,params);
|
Dawn@4
|
254
|
Dawn@4
|
255 % and call fun
|
Dawn@4
|
256 fval = feval(params.fun,xtrans,params.args{:});
|
Dawn@4
|
257
|
Dawn@4
|
258 end % sub function intrafun end
|
Dawn@4
|
259
|
Dawn@4
|
260 % ======================================
|
Dawn@4
|
261 function xtrans = xtransform(x,params)
|
Dawn@4
|
262 % converts unconstrained variables into their original domains
|
Dawn@4
|
263
|
Dawn@4
|
264 xtrans = zeros(1,params.n);
|
Dawn@4
|
265 % k allows some variables to be fixed, thus dropped from the
|
Dawn@4
|
266 % optimization.
|
Dawn@4
|
267 k=1;
|
Dawn@4
|
268 for i = 1:params.n
|
Dawn@4
|
269 switch params.BoundClass(i)
|
Dawn@4
|
270 case 1
|
Dawn@4
|
271 % lower bound only
|
Dawn@4
|
272 xtrans(i) = params.LB(i) + x(k).^2;
|
Dawn@4
|
273
|
Dawn@4
|
274 k=k+1;
|
Dawn@4
|
275 case 2
|
Dawn@4
|
276 % upper bound only
|
Dawn@4
|
277 xtrans(i) = params.UB(i) - x(k).^2;
|
Dawn@4
|
278
|
Dawn@4
|
279 k=k+1;
|
Dawn@4
|
280 case 3
|
Dawn@4
|
281 % lower and upper bounds
|
Dawn@4
|
282 xtrans(i) = (sin(x(k))+1)/2;
|
Dawn@4
|
283 xtrans(i) = xtrans(i)*(params.UB(i) - params.LB(i)) + params.LB(i);
|
Dawn@4
|
284 % just in case of any floating point problems
|
Dawn@4
|
285 xtrans(i) = max(params.LB(i),min(params.UB(i),xtrans(i)));
|
Dawn@4
|
286
|
Dawn@4
|
287 k=k+1;
|
Dawn@4
|
288 case 4
|
Dawn@4
|
289 % fixed variable, bounds are equal, set it at either bound
|
Dawn@4
|
290 xtrans(i) = params.LB(i);
|
Dawn@4
|
291 case 0
|
Dawn@4
|
292 % unconstrained variable.
|
Dawn@4
|
293 xtrans(i) = x(k);
|
Dawn@4
|
294
|
Dawn@4
|
295 k=k+1;
|
Dawn@4
|
296 end
|
Dawn@4
|
297 end
|
Dawn@4
|
298
|
Dawn@4
|
299 end % sub function xtransform end
|
Dawn@4
|
300
|
Dawn@4
|
301
|
Dawn@4
|
302
|
Dawn@4
|
303
|
Dawn@4
|
304
|