ptos_a.m

来自「FISMAT accommodates different arithmetic」· M 代码 · 共 234 行

echo off
% ptos_a.m
%
% This script runs in conjunction with the model file ptos_s.m and the 
% controller file ptos_c.m
% This script runs the simulation from workspace and analysis the results.
% 
% The block diagram is taken from Franklin, Powell, Workmann :
% 'Digital Control of Dynamic Systems'
% Addison Wesley 2nd Ed. 1990
%
% FSTB - Fuzzy Systems Toolbox for MATLAB
% Copyright (c) 1993-1996 by Olaf Wolkenhauer
% Control Systems Centre at UMIST
% Manchester M60 1QD, UK
%
% 21-May-1994

tstart=0;
tfinal=1.5;
tol=1e-3;
minstep=0.001;
maxstep=0.01;
clc
dec=input('Sequence of runs ? y/n [n]: ','s');
if isempty(dec), dec='n'; end;
if dec=='y',
  dec2=input('Saving data in file(s) ? (overwrites previous files) y/n [n]: ','s');
  if isempty(dec2)|dec2=='n'; dec2='n'; samples=[]; end;
  % final values for the reference step signal:
  ir=[1 3 4 5 6 8 9 10 -1 -3 -4 -5 -6 -8 -9 -10]; 
  % number of initial steps:
  nu_of_ir=length(ir);    
  dec3=input('Using initial steps from script file ? y/n [n]: ','s');
  if isempty(dec3), dec3='n'; end;
  if dec3=='n',
    disp('Type in a vector with initial steps like [1 2 3 4] : ')
    ir=input('');
    nu_of_ir=length(ir);
  end;
  REF=4;
  % Preparing the controller:
  echo on
  Ts=0.1
  ESET=[-10 -3 0 3 10] 
  OSET=[-45 -20 0 20 45]
  echo off
  global Ts ESET OSET REF
  dec3=input(' Changes in the controller characteristic ? y/n [n] : ','s');
  if dec3=='y',
    disp('type return to continue..');
    keyboard
  end;
  dec3=input(' Plot of the controller characteristic function ? y/n [n] : ','s');
  if dec3=='y',
    disp('Press any key after plot...')
    ctr_func('ptos_c',Ts,[min(ESET) max(ESET)]);
    title('characteristic function of the fuzzy controller');
    xlabel('error e');ylabel('u');pause
  end;
  lcount=zeros(nu_of_ir,1); % counter for the number of samples per run
  sco=0;                    % simulation counter.
  if dec2=='y',
    [infid,message]=fopen('ptos.inf','w');
    if infid==-1, disp(['Error: ',message]); end;
    fprintf(infid,'%g\n',tol);
    fprintf(infid,'%g\n',tfinal);
    fprintf(infid,'%g\n',minstep);
    fprintf(infid,'%g\n',maxstep);
    [datfid,message]=fopen('ptos.dat','w');
    if datfid==-1, disp(['Error: ',message]); end;
  end;
  clf
  axis([0 1.5 -15 15]); title('PTOS');
  hold on 
  for i=1:nu_of_ir,
    REF=ir(i);
    x0=[];
    [t,x,y]=rk45('ptos_s',[tstart,tfinal],x0,[tol,minstep,maxstep]);
    % Analysis of the results: y(1) = ctr.input y(2) = ctr.output
    % y(3) = system output
    if dec2=='y',
      fprintf(datfid,'%g %g %g %g\n',[t';y(:,1)';y(:,2)';y(:,3)']);
      spl=length(y);
      sco=sco+1;
      lcount(sco)=spl;
      fprintf(infid,'%g  %g\n',REF,spl);
    end;
    if dec2=='n',
      samples=[samples;[t,y]];
    end;
    spl=length(y);
    plot(t,y,'EraseMode','none'); %drawnow
    disp(['ir ',num2str(i)]);
  end;   % of loop for initial steps.
  hold off
  if dec2=='y',
    fprintf(infid,'%g  %g\n',sco,sum(lcount));
    status=fclose('all');
    if status==-1, disp('Error closing all files. '); end;
    disp('The data file is: ptos.dat ');
    disp('Information to these files are in ptos.inf')
  end;
  clear x0 nu_of_ir message status datfid 
  clear sco spl i ir an infid lcount t x y dec2 dec3 
end;                                       % of simulation sequence.

dec=input('Single run ? y/n [n]: ','s');
if dec=='y',
  % Preparing the controller:
  REF=4;
  echo on
  Ts=0.05
  ESET=[-10 -3 0 3 10] 
  OSET=[-45 -20 0 20 45]
  echo off
  global Ts ESET OSET REF
  dec3=input(' Changes in the controller characteristic ? y/n [n] : ','s');
  if dec3=='y',
    disp('type return to continue..');
    keyboard
  end;
  dec3=input(' Plot of the controller characteristic function ? y/n [n] : ','s');
  if dec3=='y',
    ctr_func('ptos_c',1,[min(ESET) max(ESET)]);
    title('characteristic function of the fuzzy controller');
    xlabel('error e');ylabel('u');pause
  end;
  disp(['current step height :',num2str(REF)]);
  REF=input('step height ? : ');
  x0=[];
  [t,x,y]=rk45('ptos_s',[tstart,tfinal],x0,[tol,minstep,maxstep]);
  % Analysis of the results: y(1) = ctr.input y(2) = ctr.output
  % y(3) = system output
  spl=length(y);
  clf,plot(t,y); title('ptos'); 
  %axis([0 2 0 10]);
  clear spl x0
end;               % of single simulation run.

dec=input('Loading sample data in workspace (from files) ? y/n [n]: ','s');
if dec=='y',
  % The file containing all sample vectors is called: ptos.dat .
  [datfid,message]=fopen('ptos.dat');
  if datfid==-1, disp(['Error: ',message]); end;
  samples=fscanf(datfid,'%g %g %g %g ',[4 inf]);
  samples=samples';
  status=fclose(datfid);
  if status==-1,
    [message]=ferror(datfid);
    disp(['Error: ',message]);
  end;
  disp('The simulations are now stored in "samples" .');
  clear message datfid status
end; % of loading sample data from files to workspace.
disp(' ');
disp(['There are ',num2str(length(samples)), ' sample vectors.']);
dec=input('Plotting the input data distribution ? y/n [n]: ','s');
if isempty(dec), dec='n'; end;
if dec=='y',
  clf
  plot(samples(:,2),samples(:,3),'.');
  title('Input data distribution'); xlabel('e'); ylabel('u');
end;

dec=input('Product space clustering ? y/n [n]: ','s');
if isempty(dec), dec='n'; end;
if dec=='y',
  nos=input(['Number of training pattern ? [',num2str(length(samples)),'] : ']);
  if isempty(nos), nos=length(samples); end;
  p=input('Number of synaptic vectors ? [40] : ');
  if isempty(p), p=40; end;
  N=input('Number of nearest synaptic vectors or "winners" ? [10] : ');
  if isempty(N), N=10; end;
  wii=input('wii (positive) ?  [1] : ');
  if isempty(wii), wii=1; end;
  wij=input('wij (negative) ? [-1] : ');
  if isempty(wij), wij=-1; end;
  c=input('c (c>0) ? [1] : ');
  if isempty(c), c=1; end;
  vectors=[samples(:,2) samples(:,3)];
  esec=[-17 -3;-3 0;0 3;3 17];
  usec=[-45 -20;-20 0;0 20; 20 45];
  dec=input('random mix of sample pattern ? y/n [y] : ','s');
  if isempty(dec), dec='y'; end;
  if dec=='y',
    rands=randperm(nos);
  end;
  [FAMbank,msum,y2,m2e,m2u,M,mdens]=ptos_psc(esec,usec,vectors(1:nos,:),N,p,wii,wij,c,rands);
  FAMbank,msum
  figure
  subplot(2,2,1),plot(y2); ylabel('y2');xlabel('training time');
  subplot(2,2,2),plot(m2e); ylabel('M(2,1) = e');
  subplot(2,2,4),plot(m2u); ylabel('M(2,2) = u'); xlabel('training time');
  subplot(2,2,3),axis([0 1 0 1]);set(gca,'Visible','off'); 
  text(0,0,['wii = ',num2str(wii)]);
  text(0,0.25,['wij = ',num2str(wij)]);
  text(0,.5,['pattern = ',num2str(nos)]);
  text(0,0.75,'ct = 0.1(1-t/pattern)');
  
  mdens=mdens(:).*100./p;
  figure
  subplot(2,2,1); bar(mdens);
  xlabel('cell number'); ylabel('% classification');
  title('synaptic vector histogram');
  
  subplot(2,2,2); 
  axis([esec(1,1) esec(length(esec),2) usec(1,1) usec(length(usec),2)]);
  hold on
  plot(vectors(1:nos,1),vectors(1:nos,2),'.');
  title('pattern distribution'); ylabel('u');
    
  subplot(2,2,3); axis([0 1 0 1]); set(gca,'Visible','off'); 
  text(0,0.5,['c = ',num2str(c)]);
  text(0,0.75,['p = ',num2str(p)]);
  text(0,0.25,[num2str(N),'  "winner"']);

  subplot(2,2,4); 
  axis([esec(1,1) esec(length(esec),2) usec(1,1) usec(length(usec),2)]);
  set(gca,'XTick',[esec(:,1);esec(length(esec),2)]); 
  set(gca,'XGrid','on');
  set(gca,'YTick',[usec(:,1);usec(length(usec),2)]); 
  set(gca,'YGrid','on');
  hold on; xlabel('e'); ylabel('u');
  title('Synaptic vectors M'); xlabel('e');ylabel('u');
  plot(M(1,:),M(2,:),'.'); 
  
  clear xsec esec usec vectors y2 m2e m2u p N wii wij c mdens M i
end;

clear dec tstart tfinal tol minstep maxstep nos 
echo on