Initial commit of Autotuner library in Matlab/Simulink

Autotuners/runSim.m

+% Helper program to run simulations and return resulting signals. This
+% functionality will be covered by dedicated hardware in any real world
+% implementation. The list of arguments are made available in the simulator
+% workspace.
+function [u1,y1,t,result]=runSim(varargin)
+
+% assign arguments in current workspace
+for n=1:numel(varargin)
+    feval(@(name,val)(assignin('caller',name,val)), ...
+        inputname(n),varargin{n})
+end
+
+% make current workspace available to the simulator
+options = simset('SrcWorkspace','current');
+warning('off','Simulink:blocks:TDelayDirectThroughAutoSet')
+sim('autotuner_siso_2012a',[],options)
+warning('on','Simulink:blocks:TDelayDirectThroughAutoSet')
+
+% extract signals
+t=simout.time;
+data=simout.signals.values;
+u1=data(:,1);
+y1=data(:,2);
+result = result.signals.values;
\ No newline at end of file

Autotuners/setParams.m

+
+% File setting the values of all parameters used throughout the code set
+params=[];
+
+% used by random number generator
+params.noiseSeed=0;
+
+% Noise measurement time
+params.noiseMeasureTime=10; %5
+
+% default value on minimum hysteresis level
+params.hyst = 0.1;
+%params.hyst = 0.05;
+
+%asymmetry of relay
+params.asym=2; % This one is used for the SISO autotuners
+params.gamma=1.5; % These two are used for the TITO
+params.gamma2=2;
+
+% nominal magnitude of uoff, values depends on range of control signal and
+% which one should be largest.
+params.uon = 1; % (control range 0-10)
+params.uoff = 1/params.asym; % (on-value larger than off-value)
+%params.uon = params.uoff*params.asym;
+
+
+% number of switches before experiment terminates
+params.switches=4; % SISO (We say 3 in papers, but the implementation counts the first one as a switch as well)
+%params.switches=5; % TITO
+
+% TITO-steps
+params.ystep=0;
+params.ustep=-0.8; % This parameter is set quite arbitrarily and might need further consideration. Only for TITO.
+
+% Other parameters 
+params.timeAfterExp = 0;
+params.expStartTime = 10; % Change to the required time if needed to start by moving to a stationary point
+params.ustart = 0;  % Nominal control signal. Set another level if your working point is not zero
+params.ustart2 = 0; % Nominal control signal. Set another level if your working point is not zero
+
+% sample time 
+%params.h = 0.1;
+params.h = 0.01; % for real quadtank
+
+% convergence limit
+params.eps = 0.01; % Only used for tau-tuner
+% Wanted intervals for amplitude adjustments
+params.low_limit=2; % Factor multiplied to the hysteresis level.
+params.dev_limit=12;
+% Steady-state limit
+params.equi = 0.0001;
+
+% save parameters to file
+save tmpdata/params.mat params
+
+
+% setup simulator
+simopt=[];
+simopt.h=params.h;
+simopt.uoff=-params.uoff;
+simopt.uon=-simopt.uoff*params.asym;
+simopt.ID=1; % To run relay have different from 0
+simopt.Tf = 1000; % This one should be changed to your wanted time or set to inf
+simopt.F = tf(1);
+

Autotuners/setup_script.m

+%% Run Tau-Tuner experiment
+% Runs an autotuning experiment with the Tau-tuner and saves the results to
+% file. Parameter choices are done in the file setParams and the wanted
+% autotuner and process model are chosen in the simulink model
+% autotuner_siso_2012a.mdl. 
+
+close all; clear all; clc
+% setup
+setParams
+params.expStartTime = 10; % It is possible to change the default parameters 
+%in the file setParams or to change them temporarily like this.
+
+% main experiment
+[simdata.u,simdata.y,simdata.t,simdata.result]=runSim(simopt,params);
+
+% Plot
+figure(1)
+subplot(211)
+plot(simdata.t,simdata.y)
+subplot(212)
+plot(simdata.t,simdata.u)
+
+% Save results
+s = tf('s');
+P.Kp = simdata.result(2); P.T = simdata.result(3); P.L = simdata.result(4);
+P.tf = P.Kp/(1+s*P.T)*exp(-s*P.L);
+C.K = simdata.result(5); C.Ti = simdata.result(6); C.Td = simdata.result(7);
+C.tf = C.K*(1+1/(s*C.Ti)+s*C.Td);
+margins = allmargin(P.tf);
+wo = margins.GMFrequency(1);
+tau_f = 1/(5*wo); % Filter constant
+Filter = 1/(1+s*tau_f+(s*tau_f)^2/2); % Second order filter for controller
+save('tau_tuner_results.mat', 'P', 'C','Filter','simdata', 'params')
+
+
+%% Run NOMAD experiment
+% Runs an autotuning experiment with the NOMAD-tuner and saves the results to
+% file. Parameter choices are done in the file setParams and the wanted
+% autotuner (NOMAD) and process model have to be chosen in the simulink model
+% autotuner_siso_2012a.mdl. Remember to setup_cvx before starting!
+
+close all; clear all; clc
+setParams
+addpath ../Scripts/SISO
+[simdata.u, simdata.y, simdata.t, result]=runSim(simopt,params);
+
+% Estimate models
+[P.models,idTime, P.simdata]=getPhat(params, simopt,simdata, result);
+
+% Design filter
+s=tf('s');
+model= P.models.best.P;
+margins = allmargin(model);
+wo = margins.GMFrequency(1);
+tau_f = 1/(5*wo); % Filter constant, same as in CCTA paper for TITO, but with w180 instead of wc.
+Filter = 1/(1+s*tau_f+(s*tau_f)^2/2); % Second order filter for controller
+Pfiltered = model*Filter;
+
+% Calculate controller parameters
+warning('off')
+Ms=1.4; Mt=1.4; p=[0,0,0]'; r=0;
+w=logspace(-2,2,1000)'; % Would like to base the frequency interval on the obtained model, but haven't got it to work satisfactorily yet.
+P0 = squeeze(freqresp(Pfiltered,w));
+p=cvxpid(P0,r,w,Ms,Mt,p);
+C.kp=p(1); C.ki=p(2); C.kd=p(3);
+s=tf('s');
+C.tf=C.kp+C.ki/s+C.kd*s;
+
+save('nomad_exp.mat', 'P', 'C','Filter', 'simdata', 'params','result')
+
+
+%% Run Multi-NOMAD experiment
+% Runs the autotuner_tito_2012a simulink model with the multi-NOMAD. Change
+% to your chosen process. Experiment parameters are changed in setParams.m
+  
+
+close all; clear all; clc
+setParams
+addpath ../Scripts/TITO
+simulate_experiment; % Simulates the relay exp and creates the file sim_data.mat
+identify_models; % idenitifies model from sim_data and saves to file id_models.mat
+plot_models; % Generates fit plot and bode plot
+warning('off')
+design_controllers; % Finds MIMO PID controller to model, saves to controller.mat
+%plot_controller_performance % Generates plots if you have a true model to
+%compare to
+
+clear all
+load tmpdata/sim_data.mat
+load tmpdata/id_models.mat
+load tmpdata/controllers.mat
+save('tito_exp.mat')
+
+%% Run controlled process
+% An example of how the obtained controller can be used in the SISO case.
+
+close all; clear all; clc
+s = tf('s');
+setParams
+simopt.Tf = 200; % In the simulink model a setpoint change is made after 100 s, that can be changed inside the step block.
+
+%load nomad_exp.mat
+load tau_tuner_results.mat
+FC = C.tf*Filter;
+
+options = simset('SrcWorkspace','current');
+warning('off','Simulink:blocks:TDelayDirectThroughAutoSet')
+sim('controlled_process_siso_2012a',[],options)
+warning('on','Simulink:blocks:TDelayDirectThroughAutoSet')
+simdata.u = simout.signals.values(:,1);
+simdata.y1 = simout.signals.values(:,2);
+simdata.t = simout.time;
+subplot(211)
+plot(simdata.t, simdata.y1,'b')
+ylim([0,10])
+subplot(212)
+plot(simdata.t, simdata.u,'r')
+ylim([0,10])
+save('controller_evaluation.mat','simdata')
+
+%% Run controlled process TITO
+% An example of how the obtained controller can be used in the TITO case.
+
+close all; clear all; clc
+s = tf('s');
+setParams
+simopt.Tf = 200; % Times and setpoint changes can of course be adjusted.
+
+load tito_exp.mat
+% C_hat is the MIMO controller including the second order filter!;
+C11 = C_hat.tf(1,1);
+C12 = C_hat.tf(1,2);
+C21 = C_hat.tf(2,1);
+C22 = C_hat.tf(2,2);
+options = simset('SrcWorkspace','current');
+warning('off','Simulink:blocks:TDelayDirectThroughAutoSet')
+sim('controlled_process_tito_2012a',[],options)
+warning('on','Simulink:blocks:TDelayDirectThroughAutoSet')
+simdata.u1 = simout.signals.values(:,1);
+simdata.y1 = simout.signals.values(:,2);
+simdata.u2 = simout.signals.values(:,3);
+simdata.y2 = simout.signals.values(:,4);
+simdata.t = simout.time;
+ymaxx = 5;
+subplot(411)
+plot(simdata.t, simdata.y1,'b')
+ylim([0,ymaxx])
+subplot(412)
+plot(simdata.t, simdata.u1,'r')
+ylim([0,ymaxx])
+subplot(413)
+plot(simdata.t, simdata.y2,'b')
+ylim([0,ymaxx])
+subplot(414)
+plot(simdata.t, simdata.u2,'r')
+ylim([0,ymaxx])
+save('controller_evaluation_tito.mat','simdata')
\ No newline at end of file