Commit 69be206f authored by Kristian Soltesz's avatar Kristian Soltesz
Browse files

Complete version of PIDopt

parent a1134b95
This repo contains Matlab code for optimization-based robust PID controller design.
PIDopt is Matlab code for optimization-based robust PID controller design. It comprises:
Examples, contributor credits and license information are found in PIDopt.pdf, and LICENCE, respectively.
* pidIE.m solves a load step integral error (IE) minimization problem with Hinf constraints on sensitivity, complementary sensitivity and noise sensitivity.
It yields either a PI or PID type controller. The convex-convave algorithm is described in [1].
Note that CVX [1] and MOSEK [2] must be installed and configured for full functionality.
* pidIAE.m instead solves a load integral absolute error (IAE) minimization problems. Constraints and output as with pidIE.m. It yields either a continuous-time PI or PID type controller. The exact gradient algorithm is described in [2].
[1] http://cvxr.com
[2] https://mosek.com
* pidfIAE.m has the same functionality as pidIAE.m, but simultaneously designs a second-order low-pass filter. The exact gradient algorithm is decribed in [3].
* optIAE.m solves the same problem as pidfIAE.m, but returns a high-order discrete time controller, approximating the optimal linear time invariant controller for the problem. The Youla (Q) design algorithm is described in [4].
* defaults.m generates default parameter values for the above deign scripts
* example.m uses a realistic design scenario to demonstrate the use of PIDopt.
* LICENSE, containing license information.
* README (this file).
The tool has been developed in Matlab R2016a (9.0.0.341360) 64-bit (maci64), but should work with other Matlab versions with no, or little, modification. Note that CVX [5] and MOSEK [6] must be installed and configured for full functionality.
[1]
[2]
[3]
[4]
[5] http://cvxr.com
[6] https://mosek.com
......@@ -49,7 +49,7 @@ switch optType
G = loopsens(P,K0);
[~,t] = step(G.PSi); % Load step response.
[y,t] = step(G.PSi,2*t(end)); % Longer load step response.
T = t(find(abs(y)>0.005*max(abs(y)),1,'first')); % Find 0.5% settling time.
T = t(find(abs(y)>0.005*max(abs(y)),1,'last')); % Find 0.5% settling time.
% Save to struct
opts.Nq = 50;
......
This diff is collapsed.
......@@ -33,54 +33,83 @@ opts.p = [0;0]; % Design PI controller
disp('Computing IAE-minimizing filtered PID design...')
opts = defaults(P,'pidfIAE');
opts.p = [p1;0;0]; % Hot start with PI controller
[K2,p2] = pidfIAE(P,Ms,Mt,Mks,opts);
[K2,p2,opts2] = pidfIAE(P,Ms,Mt,Mks,opts);
%% IAE-miniming Q design
disp('Computing IAE-minimizing Q design...')
opts = defaults(P,'optIAE',K2); % Hot start with filtered PID controller
[K3,Pd,Q] = optIAE(P,Ms,Mt,Mks,opts);
[K3,Pd,Q,opts3] = optIAE(P,Ms,Mt,Mks,opts);
%% Validate and visualize designs
% Pc = P;
% wN =
%
% for design = 3:-1:1 % Go backwards, to start with Q design
% P = Pc;
% eval(['K=K' num2str(design) ';'])
% switch design
% case 1
% c = 'b';
% case 2
% c = 'r';
% case 3
% c = 'g';
% P = Pd;
% wN = ; % Nyquist frequency
% end
% G = loopsens(P,K);
%
% figure(1)
%
% figure(2)
% bode(G.So,'color',c)
% end
%
%
% % opts.optIAE.K0 = K2; % Hot start with filtered PID controller.
% % %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
% %
% % % 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
% % K3 = optIAE(P,Ms,Mt,Mks,h,T,K0);
%
% %% Validate and visualize designs
\ No newline at end of file
%% Visualize and validate designs
Pc = P;
T = opts3.T; % Step response horizon
h = opts3.h; % Sampling period
wN = pi/h; % Nyquist frequency
w ={opts1.w(1),wN}; % Angular frequency san
for design = 1:3
P = Pc;
if design == 3
P = Pd;
else
P = Pc;
end
eval(['K=K' num2str(design) ';'])
G = loopsens(P,K);
figure(1)
step(G.PSi,T)
hold on
figure(2)
bode(G.So,w)
hold on
figure(3)
bode(G.To,w)
hold on
figure(4)
bode(G.CSo,w)
hold on
figure(5)
bode(K,w)
hold on
% Compute IAE
t = 0:h:T;
y = step(G.PSi,t);
IAE = trapz(t,abs(y));
eval(['IAE' num2str(design) '=IAE'])
end
% Customize plots
figure(1)
title('Load step response')
legend('IE-minimizing PI','IAE-minimizing filtered PID',...
'IAE-minimizing Q','Location','NorthEast' )
figure(2)
title('Sensitivity - S')
ax=findobj(gcf,'type','axes');
plot(ax(2),ax(2).XLim,Ms*[1 1],'k')
ylim(ax(2),[.1 5])
figure(3)
title('Complementary sensitivity - T')
ax=findobj(gcf,'type','axes');
plot(ax(2),ax(2).XLim,Mt*[1 1],'k')
ylim(ax(2),[.1 5])
figure(4)
title('Noise sensitivity - KS')
ax=findobj(gcf,'type','axes');
plot(ax(2),ax(2).XLim,Mks*[1 1],'k')
ylim(ax(2),[1 20])
figure(5)
title('Controller - K')
ax=findobj(gcf,'type','axes');
plot(ax(2),ax(2).XLim,Mks*[1 1],'k')
ylim(ax(2),[1 50])
\ No newline at end of file
function [K,Pd,Q,opts] = optIAE(P,Ms,Mt,Mks,opts)
% Computes approximation of IAE-optimal LTI controller subject to
% constraints on the magnitude of the sensitivity, complementary
% sensitivity and noise sensitivity functions.
%
% K - Approximation of the optimal LTI controller
%
% P - LTI process model.
% Ms - Maximum allowed sensitivity.
% Mt - Maximum allowed complementary sensitivity.
% Mks - Maxium allowed noise sensitivity.
% otps - Optimization options.
% Q - The Q paramter.
%
% Relies on MOSEK (https://mosek.com)
%
% Downloaded from git@gitlab.control.lth.se:kristian/PIDopt.git
% See the git repo for further documentation.
% Extract options
h = opts.h;
T = opts.T;
K0 = opts.K0;
Nq = opts.Nq;
% Sample the plant
Pd = balreal(delay2z(c2d(balreal(P),h,'foh')));
[AP,BP,CP,DP] = ssdata(Pd);
% Time horizon
Tsteps = floor(T/h); % normalize wrt samlpling period
% Optimize using CVX
cvx_clear
cvx_begin
cvx_solver mosek
cvx_precision default
cvx_quiet(true)
variable q(Nq)
% Nominal Q based on nominal controller
K0d = balreal(c2d(K0,h,'foh'));
Q0d = balreal(feedback(K0d,Pd));
if max(abs(pole(Q0d))) > 1
error('Sampled nominal design gives unstable closed loop')
end
[AQ0,BQ0,CQ0,DQ0] = ssdata(Q0d);
% Additional Q parameterization by finite pulse response
AQ1 = diag(ones(Nq-2,1),1); % Q = K/(1+PK)
BQ1 = q(2:Nq,1); % The other Nq-1 values of pulse response
CQ1 = [1 zeros(1,Nq-2)];
DQ1 = q(1); % The first (k=0) value of pulse response
[AQ,BQ,CQ,DQ] = parallell(AQ1,BQ1,CQ1,DQ1,AQ0,BQ0,CQ0,0*DQ0);
% Calculate closed-loop state-space systems
[AT,BT,CT,DT] = series(AQ,BQ,CQ,DQ,AP,BP,CP,DP); % = T = P*Q
[AS,BS,CS,DS] = parallell(0,0,0,1,AT,-BT,CT,-DT); % = S = 1-T
[APS,BPS,CPS,DPS] = series(AS,BS,CS,DS,AP,BP,CP,DP); % = P*S
if max(abs(eig(APS))) >= 1
error('P*S not stable')
end
nS = size(AS,1);
variable PS(nS,nS) symmetric
nT = size(AT,1);
variable PT(nT,nT) symmetric
if Mks < Inf
nQ = size(AQ,1);
variable PQ(nQ,nQ) symmetric
end
% Calculate pulse response d -> y
clear y
y(1) = DPS;
for k = 2:Tsteps
y(k) = CPS * APS^(k-2) * BPS;
end
y = cumsum(y); % Step response d -> y
minimize norm(y,1) % IAE
subject to:
% Constraint on S = 1/(1+PK)
XS = [PS AS*PS BS zeros(nS,1);
PS'*AS' PS zeros(nS,1) PS*CS';
BS' zeros(1,nS) 1 DS';
zeros(1,nS) CS*PS' DS Ms^2];
XS == semidefinite(2*nS+2)
% Constraint on T = PK/(1+PK)
XT = [PT AT*PT BT zeros(nT,1);
PT'*AT' PT zeros(nT,1) PT*CT';
BT' zeros(1,nT) 1 DT';
zeros(1,nT) CT*PT' DT Mt^2];
XT == semidefinite(2*nT+2)
% Constraint on Q = K/(1+PK)
if Mks < Inf
XQ = [PQ AQ*PQ BQ zeros(nQ,1);
PQ'*AQ' PQ zeros(nQ,1) PQ*CQ';
BQ' zeros(1,nQ) 1 DQ';
zeros(1,nQ) CQ*PQ' DQ Mks^2];
XQ == semidefinite(2*nQ+2)
end
% Add constraint dcgain(Q) = 1/dcgain(P);
invdcg = 1/dcgain(P);
CQ/(eye(size(AQ))-AQ)*BQ+DQ == invdcg
cvx_end
% Check for success
if strcmp(cvx_status,'Failed') || strcmp(cvx_status,...
'Inaccurate/Infeasible')
error('Optimization failed')
end
% Formulate the optimal controller
Q = ss(AQ,BQ,CQ,DQ,h);
K = feedback(Q,-Pd);
end
function [a,b,c,d] = parallell(a1,b1,c1,d1,a2,b2,c2,d2)
% Parallell interconnection of ss systems
[msg,a1,b1,c1,d1]=abcdchk(a1,b1,c1,d1); error(msg);
[msg,a2,b2,c2,d2]=abcdchk(a2,b2,c2,d2); error(msg);
[ny1,nu1] = size(d1);
[ny2,nu2] = size(d2);
% State space systems w/o selection vectors
inputs1 = [1:nu1]; outputs1 = [1:ny1];
inputs2 = [1:nu2]+nu1; outputs2 = [1:ny2]+ny1;
% Check sizes
if (length(inputs1)~=length(inputs2))
error('Input sizes don''t match.')
elseif (length(outputs1)~=length(outputs2))
error('Output sizes don''t match.')
end
% Parallel Connection
a = blkdiag(a1,a2);
b = blkdiag(b1,b2);
c = blkdiag(c1,c2);
d = blkdiag(d1,d2);
% Connect inputs
if ~isempty(b), b(:,inputs1)=b(:,inputs1)+b(:,inputs2); end
if ~isempty(d), d(:,inputs1)=d(:,inputs1)+d(:,inputs2); end
% Connect outputs
if ~isempty(c), c(outputs1,:)=c(outputs1,:)+c(outputs2,:); end
if ~isempty(d), d(outputs1,:)=d(outputs1,:)+d(outputs2,:); end
% Delete redundant inputs and outputs
b(:,inputs2) = []; d(:,inputs2) = [];
c(outputs2,:) = []; d(outputs2,:) = [];
end
\ No newline at end of file
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