opttrain.m

基于神经网络的控制工具箱

% ---------------------------------     OPTTRAIN     ------------------------------
%
%  Program for training a neural network controller as a "detuned" inverse
%  model of a process by using a method which is related to specialized training.
%
%  The neural network controller is trained to minimize the criterion
%
%                                                  
%                   ---                   2                     2
%      J(theta) =   >   [ (ym(t+1)-y(t+1))   +  rho*(u(t|theta))  ]
%                   ---                               
%                    t                             
%
%                    Bm(q)
%      where ym(t) = ----- r(t)
%                    Am(q)
%
%  The training algorithm is a recursive Gauss-Newton algorithm. Three different
%  versions have been implemented: Exponential forgetting, constant trace, and the
%  exponential forgetting and resetting algorithm (EFRA).
%
%  The user must provide a neural network model of the process and an initial
%  neural network controller. This can be created by choosing the weights at
%  random or by training a network as the inverse of the process first.
%
%  All parameters associated with the training procedure are set in the
%  file 'optrinit.m'. By choosing the penalty factor on the controls (rho) to zero
%  the function will perform exactly like 'special2' which generates an inverse
%  model by specialized training. 
%
%  (c) Magnus Norgaard IAU.
%  LastEditDate: Oct. 28, 1997. 


%----------------------------------------------------------------------------------
%-------------------         >>>  INITIALIZATIONS  <<<        ---------------------
%----------------------------------------------------------------------------------

%>>>>>>>>>>>>>>>>>>>>>>      READ VARIABLES FROM FILE       <<<<<<<<<<<<<<<<<<<<<<<
clear plot_a plot_b
global ugl
optrinit;                                % Run user specified initializations
eval(['load ' nnctrl]);                  % Load initial controller network


% >>>>>>>>>>>>>>>>>>>>>>>>   DETERMINE REGRESSOR STRUCTURE   <<<<<<<<<<<<<<<<<<<<<<   
na        = NN(1);                       % # of past y's to be used in TDL
nb        = NN(2);                       % # of past u's to be used in TDL
nk        = NN(3);                       % Time delay in system
nab       = na+sum(nb);                  % Number of "signal inputs" to each net
outputs   = 1;                           % # of outputs
inputs    = nab;                         % # of inputs
phi       = zeros(inputs+1,1);           % Initialize regressor vector


% >>>>>>>>>>>>>>>>>    DETERMINE STRUCTURE OF FORWARD MODEL    <<<<<<<<<<<<<<<<<<<<
eval(['load ' nnforw]);                  % Load forward neural model of system
hiddenf   = length(NetDeff(1,:));        % Number of hidden neurons
L_hiddenf = find(NetDeff(1,:)=='L')';    % Location of linear hidden neurons
H_hiddenf = find(NetDeff(1,:)=='H')';    % Location of tanh hidden neurons
L_outputf = find(NetDeff(2,:)=='L')';    % Location of linear output neurons
H_outputf = find(NetDeff(2,:)=='H')';    % Location of tanh output neurons
y1f       = [zeros(hiddenf,1);1];        % Hidden layer outputs
yhat      = zeros(outputs,1);            % Network output


% >>>>>>>>>>>>>>>    DETERMINE STRUCTURE OF CONTROLLER NETWORK    <<<<<<<<<<<<<<<<
hiddenc   = length(NetDefc(1,:));        % Number of hidden neurons
L_hiddenc = find(NetDefc(1,:)=='L')';    % Location of linear hidden neurons
H_hiddenc = find(NetDefc(1,:)=='H')';    % Location of tanh hidden neurons
L_outputc = find(NetDefc(2,:)=='L')';    % Location of linear output neurons
H_outputc = find(NetDefc(2,:)=='H')';    % Location of tanh output neurons
y1c       = [zeros(hiddenc,1);1];        % Hidden layer outputs
d21       = W2f(1:hiddenf);              % Derivative if linear output
d10       = W1f(:,na+1);                 % Derivative if linear hidden units
d20       = 0;                           % Derivative of output w.r.t. control


%>>>>>>>>>>>>>>>>>    CALCULATE REFERENCE SIGNAL & FILTER IT     <<<<<<<<<<<<<<<<<<
if strcmp(refty,'siggener'),
  ref = zeros(samples+1,1);
  for ii = 1:samples,
    ref(ii) = siggener(Ts*(ii-1),sq_amp,sq_freq,sin_amp,sin_freq,dc,sqrt(Nvar));
  end
else
  eval(['ref = ' refty ';']);
  ref=ref(:);
  i=length(ref);
  if i>samples+1,
    ref=ref(1:samples+1);
  else
    ref=[ref;ref(i)*ones(samples+1-i,1)];
  end
end
ym = filter(Bm,Am,ref);               % Filter the reference
ym(samples+1) = ym(1);                % Necessary because the reference is repeated
ref(samples+1) = ref(1);


%>>>>>>>>>>>>>>>>>>>>>>>        INITIALIZE VARIABLES        <<<<<<<<<<<<<<<<<<<<<<<
% Initialization of vectors containing past signals
maxlength = 5;
y_old     = zeros(maxlength,1);
u_old     = zeros(maxlength,1);

% Variables associated with the weight update algorithm
SSE = 0;                                % Sum of squared error in current epoch
J   = 0;                                % Value of cost function in current epoch
epochs = 0;                             % Epoch counter
first = max(na,nb+nk-1)+10;             % Update weights when iteration>first
index = (hiddenc+1) + 1 + [0:hiddenc-1]*(inputs+1); % A useful vector!
parameters1= hiddenc*(inputs+1);        % # of input-to-hidden weights
parameters2= (hiddenc+1);               % # of hidden-to-output weights
parameters = parameters1 + parameters2; % Total # of weights
PSI        = zeros(parameters,1);       % Deriv. of each output w.r.t. each weight
                                        % Parametervector containing all weights
theta = [reshape(W2c',parameters2,1) ; reshape(W1c',parameters1,1)];
index3= 1:(parameters+1):(parameters^2);% Yet another useful vector
if strcmp(method,'ff'),                 % Forgetting factor method
  mflag     = 1;                        % Method flag
  lambda    = trparms(1);               % Forgetting factor
  p0        = trparms(2);               % Diagonal element of covariance matrix
  P         = p0 * eye(parameters);     % Initialize covariance matrix
elseif strcmp(method,'ct'),             % Constant trace method
  mflag     = 2;                        % Method flag
  lambda    = trparms(1);               % Forgetting factor
  alpha_max = trparms(2);               % Max. eigenvalue
  alpha_min = trparms(3);               % Min. eigenvalue
  P      = alpha_max * eye(parameters); % Initialize covariance matrix
  Pbar      = P;
elseif strcmp(method,'efra'),           % EFRA method
  mflag     = 3;                        % Method flag
  alpha     = trparms(1);               % EFRA parameters
  beta      = trparms(2);
  delta     = trparms(3);
  lambda    = trparms(4);
  gamma     = (1-lambda)/lambda;
  maxeig = gamma/(2*delta)*(1+sqrt(1+4*beta*delta/(gamma*gamma)));% Max. eigenvalue
  P      = maxeig * eye(parameters);    % Initialize covariance matrix
  betaI     = beta*eye(parameters);     % Useful diagonal matrix
end

% Initialization of Simulink system
if strcmp(simul,'simulink')
  simoptions = simset('Solver',integrator,'MaxRows',0); % Set integrator opt.
  eval(['[sizes,x0] = ' sim_model '([],[],[],0);']);    % Get initial states
end

% Initializations of vectors used for storing old data
ref_data    = [ref(1:samples)];
ym_data     = [ym(1:samples)];
u_data      = zeros(samples,1);
y_data      = zeros(samples,1);
yhat_data   = zeros(samples,1);

% Miscellanous initializations
maxiter = maxiter*samples;              % Number of iterations
u   = 0;
y   = 0;
t   = -Ts;
i   = 0;                                % Iteration in current epoch counter
fighandle=progress;

%----------------------------------------------------------------------------------
%---------------------         >>>   MAIN LOOP   <<<           --------------------
%----------------------------------------------------------------------------------
for iter=1:maxiter,
  i = i+1;
  t = t + Ts;


  %>>>>>>>>>>>>>>  PREDICT OUTPUT OF SYSTEM USING THE FORWARD MODEL   <<<<<<<<<<<<<
  phi = [y_old(1:na);u_old(1:nb);1];
  h1f = W1f*phi;
  y1f(H_hiddenf)  = pmntanh(h1f(H_hiddenf));
  y1f(L_hiddenf)  = h1f(L_hiddenf);    
  h2f = W2f*y1f;
  yhat(H_outputf) = pmntanh(h2f(H_outputf));
  yhat(L_outputf) = h2f(L_outputf);


  %>>>>>>>>>>>>>>>>>>>>  READ OUTPUT FROM THE PHYSICAL SYSTEM   <<<<<<<<<<<<<<<<<<<
  if strcmp(simul,'simulink')
    utmp=[t-Ts,u_old(1);t,u_old(1)];
    simoptions.InitialState=x0;
    [time,x0,y] = sim(sim_model,[t-Ts t],simoptions,utmp);
    x0 = x0(size(x0,1),:)';
    y  = y(size(y,1),:)';
  elseif strcmp(simul,'matlab')
    ugl = u_old(1);
    [time,x] = ode45(mat_model,[t-Ts t],x0);