Commit 09e6f04e authored by Kristian Soltesz's avatar Kristian Soltesz
Browse files

Updated core code and added Q design. example.m is not complete yet and...

Updated core code and added Q design. example.m is not complete yet and optIAE.m needs to be validated and integrated with defauults.m
parent f5d2fb83
This repo contains Matlab code for optimization-based robust PID controller design.
Examples, contributor credits and license information are found in PIDopt.pdf, and LICENCE, respectively.
Note that CVX [1] and MOSEK [2] must be installed and configured for full functionality.
[1] http://cvxr.com
[2] https://mosek.com
function w = defaultGrid(P)
% Get default angular frequency grid based on process dynamics.
%
% P - LTI process model.
%
% Downloaded from git@gitlab.control.lth.se:kristian/PIDopt.git
% See the git repo for further documentation.
N = 5e2; % Number of grid points
[~,~,w]=bode(P); % Grequenies Matlab considers interesting
wMax = w(end);
wMin = min(w(1),wMax/100); % ensure to include low frequencies
w=logspace(log10(wMin),log10(wMax),N)';
\ No newline at end of file
function opts = defaults(P)
% Returns default configuration options for the design scripts of PIDopt.
%
% P - Process dynamics. Used to determine angular frequency grid.
%
% Downloaded from git@gitlab.control.lth.se:kristian/PIDopt.git
% See the git repo for further documentation.
% angular frequency grid
Nw = 1e3; % number of grid points
[~,~,w] = bode(P); % generate relevant range
w = logspace(log10(w(1))-2,log10(w(end))+2,Nw)'; % extend range
% pidIE
pidIE.w = w; % Angular frequency grid for cont straint evaluation.
pidIE.p = [0;0;0]; % Vector of controller initial controller parameters
% [kp;ki;kd]. Length must be 2 (PI) or 3 (PID).
% Controller must stabilize plant.
pidIE.maxIter = 100; % Number of convex-concave iterations
pidIE.stopTol = 1e-3; % Relative objective decrease stop criterium
% pidIAE
pidIAE.w = w; % Angular frequency grid for constraint evaluation.
pidIAE.p = [0;0;0]; % Vector of initial controller parameters [kp;ki;kd].
% Length must be 2 (PI) or 3 (PID).
% Controller must stabilize plant.
pidIAE.fminOpts = optimset('algorithm','active-set','GradConstr','on',...
'GradObj','on','Display','off','maxIter',10); % fminCon options
% pidfIAE
pidfIAE.w = w; % Angular frequency grid for constraint evaluation.
pidfIAE.pC = [0;0;0]; % Vector of initial controller parameters [kp;ki;kd].
% Controller must stabilize plant.
pidfIAE.pF = 0; % Initial filter time constant.
pidfIAE.fminOpts = optimset('algorithm','active-set','GradConstr','on',...
'GradObj','on','Display','off','maxIter',50); % fmincon options
% optIAE
optIAE = []; % FIXME: implement
opts.pidIE = pidIE;
opts.pidIAE = pidIAE;
opts.pidfIAE = pidfIAE;
opts.optIAE = optIAE;
\ No newline at end of file
% This example demonstrates how to invoke the scripts for obtaining
% IE-optimal and IAE optimal (with or without LP filter) PID controllers.
% This example demonstrates how to perform a design using PIDopt.
% See defaults.m for a complete list of configurable options. (It might for
% example be that your design problem requires a custom angular frequency
% grid.)
%
% Downloaded from git@gitlab.control.lth.se:kristian/PIDopt.git
% See the git repo for further documentation.
%% Setup
clear all
clc
clear all; close all;clc % start from blank sheet
% Process model (continuous time LTI SISO system)
% Process model (continuous time LTI SISO system).
s = zpk('s');
P = 1/(s^2+2*.5*s+1)*exp(-0.5*s);
T = 4; % Time constant
L = 1; % Delay
K = 1; % Gain
P = K/(s*T+1)*exp(-s*L); % priocess dynamics
% Robustness consraints
Ms = 1.8; % Maximum allowed sensitivity function magnitude
Mt = Ms; % Maximum allowed complementary sensitivity function magnitude
Mks = 2; % Maximum allowed magnitude of transfer function from process
Ms = 1.5; % Maximum allowed sensitivity function magnitude
Mt = Ms; % Maximum allowed complementary sensitivity function magnitude
Mks = 10; % Maximum allowed magnitude of transfer function from process
% output to control signal, sometimes referred to as noise
% sensitivity.
%% Design
p0 = [0 0 0]'; % Parameters of any PID controller, which stabilizes P.
% See pidIE.m for details. Use p = [0 0] to synthesize PI
% controller
w = defaultGrid(P); % default angular frequency grid
[K1,p1] = pidIE(P,Ms,Mt,Mks,p0,w) % IE-optial PID without filter
p0 = p1; % hot start
[K2,p2] = pidIAE(P,Ms,Mt,Mks,p1,w) % IAE-optimal PID withot filter
p0 = p2; % hot start
[K3,p3] = pidfIAE(P,Ms,Mt,Mks,p2,w) % IAE-optimal PID with filter
%% Evaluation
% Closed-loop transfer functions
L = @(K)series(P,K);
S = @(K)feedback(1,L(K));
T = @(K)feedback(L(K),1);
KS = @(K)feedback(K,P);
PS = @(K)feedback(P,K);
wMin = w(1);
wMax = w(end);
close all
% Controller dynamics
figure(1)
bode(K1,{wMin,wMax},'b',K2,'r',K3,'g')
hold on
title('Controller')
legend('IE-optimal','IAE-optimal','IAE-optimal with filter','Location',...
'NorthWest')
% Bode magnitudes (to verify robustness constraints)
figure(2)
subplot(311)
bodemag(S(K1),{wMin,wMax},'b',S(K2),'r',S(K3),'g')
hold on
ylim([1e-1 2*Ms])
set(gca,'ytick',[1 Ms])
plot([wMin wMax],Ms*[1 1],'k')
title('Sensitivity magnitude')
legend('IE-optimal','IAE-optimal','IAE-optimal with filter','Constraint',...
'Location','NorthWest')
subplot(312)
bodemag(T(K1),{wMin,wMax},'b',T(K2),'r',T(K3),'g')
hold on
ylim([1e-1 2*Mt])
set(gca,'ytick',[1 Mt])
plot([wMin wMax],Mt*[1 1],'k')
title('Complementary sensitivity')
subplot(313)
bodemag(KS(K1),{wMin,wMax},'b',KS(K2),'r',KS(K3),'g')
hold on
ylim([1e-1 2*Mks])
set(gca,'ytick',[1 Mks])
plot([wMin wMax],Mks*[1 1],'k')
title('Noise sensitivity')
% Load disturbance response (to compute IE and IAE)
figure(3)
step(PS(K1),'b',PS(K2),'r',PS(K3),'g',P,'k')
hold on
title('Load step response')
legend('IE-optimal','IAE-optimal','IAE-optimal with filter','Open-loop',...
'Location','East')
\ No newline at end of file
opts = defaults(P); % Default options
%% IE-minimizing PI design
opts.pidIE.p = [0;0]; % design PI controller
[K1,p1,opts1] = pidIE(P,Ms,Mt,Mks,opts.pidIE);
%% IAE-minimizing filtered PID design
opts.pidfIAE.p = [p1;0;0]; % Hot start with PI controller.
[K2,p2] = pidfIAE(P,Ms,Mt,Mks,opts.pidfIAE); % IAE-optimal PID with filter
%% IAE-miniming Q design
K0 = K2;
h = 0.1/bandwidth(feedback(K0*P,1));
if totaldelay(P)/h > 15 % avoid excessive number of states for delay-dominated plants
h = totaldelay(P)/15;
end
h
% Determine time horizon for optimization
loops = loopsens(P,K0); % GoF
[y,t] = step(series(P,loops.So)); % load step response
T0 = t(end);
[y,t] = step(series(P,loops.So),2*T0); % longer load step response
iy = find(abs(y)>0.005*max(abs(y))); % Find 0.5% settling time
T = t(iy(end))
% Optimize
'Computing IAE optimal LTI controller'
K3 = optIAE(P,Ms,Mt,Mks,h,T,K0);
%% Validate and visualize designs
\ No newline at end of file
% The relative dead time tau = L/(L+T) and arrest time Tar = T+L are
% defined for first order time SISO LTI systems K/(s*T+1)*exp(-s*L).
%
% This code computes generalizations of tau and Tar for arbitrary SISO LTI
% systems, which can be either discrete or continuous.
%
% P - LTI SISO system.
%
% Downloaded from git@gitlab.control.lth.se:kristian/PIDopt.git
function [tau,Tar,L] = getTauTar(P)
[~,t] = step(P);
t = linspace(0,t(end),1e5);
[y,t] = step(P,t);
ydot = impulse(P,t);
[~,i] = max(ydot);
tmax = t(i);
ymax = y(i);
L = tmax-ymax/ydot(i);
if dcgain(P) == inf
tau = 0;
Tar = inf;
return
end
P = P/dcgain(P);
if P.Ts > 0
P = d2c(P,'zoh')
end
Tar = sum(-1./pole(P))-sum(-1./zero(P))+totaldelay(P);
T = Tar-L;
tau = L/Tar;
\ No newline at end of file
function [K,p,w] = pidIAE(P,Ms,Mt,Mks,p,w)
function [K,p,opts] = pidIAE(P,Ms,Mt,Mks,opts)
% Computes PID controller that maximizes the integral gain subject to
% constraints on the magnitude of the sensitivity, complementary
% sensitivity and noise sensitivity fu nctions.
% The length of p decides if the controller is of PI or PID type.
% sensitivity and noise sensitivity functions.
%
% K(s) = kp + ki/s + kd*s
%
% P - LTI process model.
% Ms - Maximum allowed sensitivity.
% Mt - Maximum allowed complementary sensitivity.
% Mks - Maxium allowed noise sensitivity.
% p - Vector of initial controller parameters. Length must be 2 (PI)
% or 3 (PID). Controller must stabilize plant.
% w - angular freqency grid. Use [] for default grid.
% P - LTI process model.
% Ms - Maximum allowed sensitivity.
% Mt - Maximum allowed complementary sensitivity.
% Mks - Maxium allowed noise sensitivity.
% otps - Optimization options. Use [] to generate defaults.
% p - Parameter vector: [kp;ki;kd] for PID controller.
%
% Downloaded from git@gitlab.control.lth.se:kristian/PIDopt.git
% See the git repo for further documentation.
% Generate default options if not provided as argument
if isempty(opts)
opts = defaults(P);
opts = opts.pidIAE;
end
% Extract optsions
w = opts.w;
p = opts.p;
fminOpts = opts.fminOpts;
% regularize
if p(1) ~= 0
k0 = p(1);
......@@ -26,11 +35,6 @@ else
k0 = 1;
end
% Frequency grid (heuristic) - change as needed
if isempty(w)
w = defaultGrid(P);
end
% Helper functions
mdot = @(x,y)(bsxfun(@times,x,y)); % (row vector)-matrix column-wise dot
fr = @(x)permute(freqresp(x,w),[3 1 2]); % frequency response over grid
......@@ -48,9 +52,7 @@ Pw = fr(P);
K_pw = fr(K_p);
% Optimization
cfg = optimset('algorithm','active-set','GradConstr','on','GradObj','on',...
'Display','off');
[p,IAE] = fmincon(@(p)J(p),p,[],[],[],[],eps+0*p,[],@(p)cNL(p),cfg);
[p,IAE] = fmincon(@(p)J(p),p,[],[],[],[],eps+0*p,[],@(p)cNL(p),fminOpts);
p = p*k0;
K = K(p);
......
function [K,p,w] = pidIE(P,Ms,Mt,Mks,p,w)
function [K,p,opts] = pidIE(P,Ms,Mt,Mks,opts)
% Computes PID controller that maximizes the integral gain subject to
% constraints on the magnitude of the sensitivity, complementary
% sensitivity and noise sensitivity functions.
% The length of p decides if the controller is of PI or PID type.
%
% K(s) = kp + ki/s + kd*s
%
% P - LTI process model.
% Ms - Maximum allowed sensitivity.
% Mt - Maximum allowed complementary sensitivity.
% Mks - Maxium allowed noise sensitivity.
% p - Vector of controller parameters. Length must be 2 (PI)
% or 3 (PID). Controller must stabilize plant.
% w - angular freqency grid. Use [] for default grid.
% P - LTI process model.
% Ms - Maximum allowed sensitivity.
% Mt - Maximum allowed complementary sensitivity.
% Mks - Maxium allowed noise sensitivity.
% otps - Optimization options. Use [] to generate defaults.
% p - Parameter vector: [kp;ki;kd] for PID controller.
%
% Relies on CVX (http://cvxr.com)
%
% Downloaded from git@gitlab.control.lth.se:kristian/PIDopt.git
% See the git repo for further documentation.
% Frequency grid (heuristic) - change as needed
if isempty(w)
w = defaultGrid(P);
% Generate default options if not provided as argument
if isempty(opts)
opts = defaults(P);
opts = opts.pidIE;
end
% Extract optsions
w = opts.w;
p = opts.p;
maxIter = opts.maxIter; % 100
stopTol = opts.stopTol; % 1e-3
% Frequency response
P = squeeze(freqresp(balreal(P),w));
......@@ -38,9 +45,9 @@ K1=K1(:,1:n);
K=K1*p;
% Optimization routine
prev_obj = -inf; % Used for the stopping criterion
max_iters = 10; % Max number of iterations (heuristic)
for iter = 1:max_iters
prevObj = -inf; % Used for the stopping criterion
maxIter = 100; % Max number of iterations (heuristic)
for iter = 1:maxIter
Lc = P.*K; % Loop transfer funtion
cvx_begin
cvx_quiet(true)
......@@ -59,15 +66,15 @@ for iter = 1:max_iters
abs(K) - Mks.*real(conj(1+Lc)./(abs(1+Lc)).*(1+L)) <= 0; % Noise sensitivity constraint
end
cvx_end
% Stopping criterion (heuristic) - change as needed
if (p(2)-prev_obj) < 1e-3*prev_obj
if (p(2)-prevObj) < stopTol*prevObj
break;
end
prev_obj = p(2);
prevObj = p(2);
end
s = zpk('s');
K1 = [1 1/s s];
K1 = K1(1:n);
K = K1*p;
K = K1*p;
\ No newline at end of file
function [K,p,w]=pidfIAE(P,Ms,Mt,Mks,pC,w)
function [K,p,opts] = pidfIAE(P,Ms,Mt,Mks,opts)
% Computes PID parameters, that minimize load step IAE, and filter for
% noise attenuation.
%
% K(s) = F(s) * C(s) = (kp + ki/s + kd*s) * 1/((s*T)^2+2*z*T*s+1)
%
% P - LTI process model.
% Ms - Maximum allowed sensitivity.
% Mt - Maximum allowed complementary sensitivity.
% Mks - Maxium allowed noise sensitivity. This is only meaningful if
% there is high-frequency roll-off (through e.g. a filter).
% pC - Vector of initial controller parameters [kp ki kd].
% Length must be 2 (PI)
% or 3 (PID). Controller must stabilize plant.
% P - LTI process model.
% Ms - Maximum allowed sensitivity.
% Mt - Maximum allowed complementary sensitivity.
% Mks - Maxium allowed noise sensitivity. This is only meaningful if
% there is high-frequency roll-off (through e.g. a filter).
% otps - Optimization options. Use [] to generate defaults.
% p - Parameter vector: [kp;ki;kd;T]
%
% Downloaded from git@gitlab.control.lth.se:kristian/PIDopt.git
% See the git repo for further documentation.
% Omit designing filter if Mks=Inf
if Mks == Inf
opts.p = opts.p(1:end-1);
[K,p,opts] = pidIAE(P,Ms,Mt,Mks,opts);
T = 0;
p = [p;zeros(4-numel(p))];
return
end
% Generate default options if not provided as argument
if isempty(opts)
opts = defaults(P);
opts = opts.pidIAE;
end
% Extract optsions
w = opts.w;
pC = opts.pC;
pF = opts.pF;
fminOpts = opts.fminOpts;
% regularize
if pC(1) ~= 0
c0 = pC(1);
......@@ -25,14 +45,9 @@ else
c0 = 1;
end
% Frequency grid (heuristic) - change as needed
if isempty(w)
w = defaultGrid(P);
end
% Ensure LP filter breakdown within grid
Tmin = 1/w(end);
pF = Tmin;
pF = max(pF,Tmin);
% Organize parameters parameters
p = [pC(:);pF(:)];
......@@ -65,9 +80,7 @@ K_pw = @(p)fr(K_p(p));
% Optimization
cfg = optimset('algorithm','active-set','GradConstr','on','GradObj',...
'on','Display','off','maxIter',10);
p = fmincon(@(p)J(p),p,[],[],[],[],[eps 0 0 Tmin],[],@(p)cNL(p),cfg);
p = fmincon(@(p)J(p),p,[],[],[],[],[eps 0 0 Tmin],[],@(p)cNL(p),fminOpts);
p(1:end-1) = p(1:end-1)*c0;
K = minreal(K(p));
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment