www.gusucode.com > nncontrol 工具箱 matlab 源码程序 > nncontrol/csrchhyb.m
function [up_delta,J,dJdu_old,dJdu,retcode,delta,tol] = csrchhyb(up,u_vec,ref,Ai,Nu,N1,N2,d,Ni,Nj,dX, ...
dJdu,J,dperf,delta,rho,dUtilde_dU,alpha,tol,Ts,min_i,max_i,Normalize,minp,maxp)
%CSRCHHYB One-dimensional minimization using a hybrid bisection-cubic search.
%
% Syntax
%
% [up_delta,J,dJdu_old,dJdu,retcode,delta,tol] = csrchhyb(up,u_vec,ref,Ai,Nu,N1,N2,d,Ni,Nj,dX, ...
% dJdu,J,dperf,delta,rho,dUtilde_dU,alpha,tol,Ts,min_i,max_i,Normalize)
%
% Description
%
% CSRCHHYB is a linear search routine. It searches in a given direction
% to locate the minimum of the performance function in that direction.
% It uses a technique which is a combination of a bisection and a
% cubic interpolation.
%
% CSRCHHYB(...) takes these inputs,
% up - Plant Inputs during the Control Horizon (Nu).
% u - Plant Inputs during the Cost Horizon (N2).
% ref - Reference input.
% Ai - Initial input delay conditions.
% Nu - Control Horizon.
% N1 - Beginning of the Control and Cost Horizons (Usually 1).
% N2 - Cost Horizon.
% d - Counter that defined initial time (Usually 1).
% Ni - Number of delayed plant inputs.
% Nj - Number of delayed plant outputs.
% dX - Search direction vector for U.
% dJdu - Derivate of the cost function respect U.
% J - Cost function value.
% dperfa - Slope of performance value at current U in direction of dX.
% delta - Initial step size.
% rho - Control weighting factor.
% dUtlde_dU - Derivate of the difference of U(t)-U(t-1) respect U.
% alpha - Search parameter.
% tol - Tolerance on search.
% Ts - Time steps.
% min_i - Minimum Input to the Plant.
% max_i - Maximum Input to the Plant.
% Normalize - Indicate if the NN has input-output normalized.
% and returns,
% up_delta - New Plant Inputs for the Control Horizon (Nu).
% J - New Cost function value.
% dJdu_old - Previous Derivate of the cost function respect U.
% dJdu - New Derivate of the cost function respect U.
% RETCODE - Return code which has three elements. The first two elements correspond to
% the number of function evaluations in the two stages of the search
% The third element is a return code. These will have different meanings
% for different search algorithms. Some may not be used in this function.
% 0 - normal; 1 - minimum step taken; 2 - maximum step taken;
% 3 - beta condition not met.
% DELTA - New initial step size. Based on the current step size.
% TOL - New tolerance on search.
%
% Parameters used for the hybrid bisection-cubic algorithm are:
% alpha - Scale factor which determines sufficient reduction in perf.
% beta - Scale factor which determines sufficiently large step size.
% bmax - Largest step size.
% scale_tol - Parameter which relates the tolerance tol to the initial step
% size delta. Usually set to 20.
% The defaults for these parameters are set in the training function which
% calls it. See TRAINCGF, TRAINCGB, TRAINCGP, TRAINBFG, TRAINOSS
%
% Algorithm
%
% CSRCHHYB locates the minimum of the performance function in
% the search direction dX, using the hybrid
% bisection-cubic interpolation algorithm described on page 50 of Scales.
% (Introduction to Non-Linear Estimation 1985)
%
% See also CSRCHBAC, CSRCHBRE, CSRCHCHA, CSRCHGOL
%
% References
%
% Scales, Introduction to Non-Linear Estimation, 1985.
% Orlando De Jesus, Martin Hagan, 1-30-00
% Copyright 1992-2002 The MathWorks, Inc.
tiu = d-N1+Ni;
upi = [1:Nu-1 Nu(ones(1,N2-d-Nu+2))]; % [1 2 ... Nu Nu ... Nu]
uvi = [tiu:N2-N1+Ni];
u = 999.9;
perfu = 999.99;
dperfu = 999.99;
% ALGORITHM PARAMETERS
scale_tol = 20;
beta = 0.9;
bmax = 26;
min_grad = 1e-6;
delta_orig=delta;
% STEP SIZE INCREASE FACTOR FOR INTERVAL LOCATION (NORMALLY 2)
scale = 2;
% INITIALIZE A AND B
a = 0;
a_old = 0;
b = delta;
perfa = J;
dperfa = dperf;
perfa_old = perfa;
dperfa_old = dperfa;
cnt1 = 0;
cnt2 = 0;
up_delta = max(min(up + b*dX,max_i),min_i); % A priori iteration
u_vec(uvi) = up_delta(upi); % Insert updated controls
% CALCLULATE PERFORMANCE FOR B
[JJ,dJJ]=calcjjdjj(u_vec,Ni,Nu,Nj,N2,Ai,Ts,ref,tiu,rho,dUtilde_dU,Normalize,minp,maxp);
perfb = JJ;
dJdub = dJJ;
dperfb = dJdub'*dX;
cnt1 = cnt1 + 1;
% INTERVAL LOCATION
% FIND INITIAL INTERVAL WHERE MINIMUM PERF OCCURS
while (perfa>perfb)&(b<bmax)
a_old=a;
perfa_old = perfa;
dperfa_old = dperfa;
perfa = perfb;
dperfa = dperfb;
a=b;
b=scale*b;
up_delta = max(min(up + b*dX,max_i),min_i); % A priori iteration
u_vec(uvi) = up_delta(upi); % Insert updated controls
% CALCLULATE PERFORMANCE FOR B
[JJ,dJJ]=calcjjdjj(u_vec,Ni,Nu,Nj,N2,Ai,Ts,ref,tiu,rho,dUtilde_dU,Normalize,minp,maxp);
perfb = JJ;
dJdub = dJJ;
dperfb = dJdub'*dX;
cnt1 = cnt1 + 1;
end
if (a == a_old)
% TAKE INITIAL BISECTION STEP IF NO MIDPOINT EXISTS
x = (a + b)/2;
X_step = x*dX;
up_delta = max(min(up + X_step,max_i),min_i); % A priori iteration
u_vec(uvi) = up_delta(upi); % Insert updated controls
% CALCLULATE PERFORMANCE FOR X
[JJ,dJJ]=calcjjdjj(u_vec,Ni,Nu,Nj,N2,Ai,Ts,ref,tiu,rho,dUtilde_dU,Normalize,minp,maxp);
perfx = JJ;
dJdux = dJJ;
dperfx = dJdux'*dX;
cnt1 = cnt1 + 1;
else
% USE ALREADY COMPUTED VALUE AS INITIAL BISECTION STEP
x = a;
perfx = perfa;
dperfx = dperfa;
a=a_old;
perfa=perfa_old;
dperfa = dperfa_old;
end
% DETERMINE THE W POINT (A OR B WITH MINIMUM FUNCTION VALUE)
if perfa>perfb
w = b;
perfw = perfb;
dperfw = dperfb;
else
w = a;
perfw = perfa;
dperfw = dperfa;
end
% DETERMINE THE OVERALL MINIMUM POINT
minperf = min([perfa perfb perfx]);
amin = a; dperfmin = dperfa;
if perfb<= minperf
amin = b; dperfmin = dperfb;
elseif perfx <= minperf
amin = x; dperfmin = dperfx;
end
% LOCATE THE MINIMUM POINT BY THE HYBRID BISECTION-CUBIC SEARCH
while ((b-a)>tol) & ((minperf > J + alpha*amin*dperf) | abs(dperfmin)>abs(beta*dperf) )
if(abs(w-x)<.02*(b-a))
bisection = 1;
else
% CUBIC INTERPOLATION
if (w > x)
aa = x; fa = perfx; ga = dperfx;
bb = w; fb = perfw; gb = dperfw;
else
bb = x; fb = perfx; gb = dperfx;
aa = w; fa = perfw; ga = dperfw;
end
ww = 3*(fa - fb)/(bb-aa) + ga + gb;
w_gagb = ww^2 - ga*gb;
if (w_gagb >= 0)
v = sqrt(w_gagb);
u_star = aa + (bb-aa)*(1 - (gb + v - ww)/(gb - ga +2*v));
if ((u_star > a)&(u_star < b))
up_delta = max(min(up + u_star*dX,max_i),min_i); % A priori iteration
u_vec(uvi) = up_delta(upi); % Insert updated controls
% CALCLULATE PERFORMANCE FOR U_STAR
[JJ,dJJ]=calcjjdjj(u_vec,Ni,Nu,Nj,N2,Ai,Ts,ref,tiu,rho,dUtilde_dU,Normalize,minp,maxp);
perfu = JJ;
dJduu = dJJ;
dperfu = dJduu'*dX;
u = u_star;
cnt2 = cnt2 + 1;
bisection = 0;
else
bisection = 1;
end
else
bisection = 1;
end
end
if (bisection == 1)
% BISECTION
if ((dperfa<0) & ((dperfx>0) | (perfx>perfa))) | ((dperfa>0) & (dperfx>0) & (perfx<perfa))
u = (a + x)/2;
else
u = (x + b)/2;
end
up_delta = max(min(up + u*dX,max_i),min_i); % A priori iteration
u_vec(uvi) = up_delta(upi); % Insert updated controls
% CALCLULATE PERFORMANCE FOR U
[JJ,dJJ]=calcjjdjj(u_vec,Ni,Nu,Nj,N2,Ai,Ts,ref,tiu,rho,dUtilde_dU,Normalize,minp,maxp);
perfu = JJ;
dJduu = dJJ;
dperfu = dJduu'*dX;
cnt2 = cnt2 + 1;
end
if ( dperfu < min_grad )
a = u; perfa = perfu; dperfa = dperfu;
b = u; perfb = perfu; dperfb = dperfu;
elseif (u>x)
a = x; perfa = perfx; dperfa = dperfx;
elseif (u<x)
b = x; perfb = perfx; dperfb = dperfx;
else
a = x; perfa = perfx; dperfa = dperfx;
b = x; perfb = perfx; dperfb = dperfx;
end
% DETERMINE THE W POINT (A OR B WITH MINIMUM FUNCTION VALUE)
if perfa>perfb
w = b;
perfw = perfb;
dperfw = dperfb;
else
w = a;
perfw = perfa;
dperfw = dperfa;
end
x = u; perfx = perfu; dperfx = dperfu;
minperf = min([perfa perfb perfx]);
amin = a; dperfmin = dperfa;
if perfb<= minperf
amin = b; dperfmin = dperfb;
elseif perfx <= minperf
amin = x; dperfmin = dperfx;
end
end
perf = minperf;
a = amin;
% COMPUTE FINAL GRADIENT
up_delta = max(min(up + a*dX,max_i),min_i); % A priori iteration
u_vec(uvi) = up_delta(upi); % Insert updated controls
% CALCLULATE PERFORMANCE FOR FINAL GRADIENT
[JJ,dJJ]=calcjjdjj(u_vec,Ni,Nu,Nj,N2,Ai,Ts,ref,tiu,rho,dUtilde_dU,Normalize,minp,maxp);
J = JJ;
dJdu_old=dJdu;
dJdu = dJJ;
% CHANGE INITIAL STEP SIZE TO PREVIOUS STEP
delta=a;
if delta < delta_orig
delta = delta_orig;
end
if tol>delta/scale_tol
tol=delta/scale_tol;
end
retcode = [cnt1 cnt2 0];