Commit 71a8178e authored by Kristian Soltesz's avatar Kristian Soltesz
Browse files

Done with core code. Tagging this version 1

parent 69be206f
......@@ -19,9 +19,17 @@ It yields either a PI or PID type controller. The convex-convave algorithm is de
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]
References:
[1] M. Hast, K. J. Astrom, B. Bernhardsson, S. Boyd. PID design by convex-concave optimization. European Control Conference. IEEE. Zurich, Switzerland. 2013.
[2] C. Grimholt, S. Skogesatd. Improved Optimization-based Design of PID Controllers Using Exact Gradients. 12th International Symposium on Process Systems Engineering and 25th European Symposium on Computer Aided Process Engineering. Elsevier. Copenhagen, Denmark. 2015.
[3] K. Soltesz, C. Grimholt, S. Skogestad. Simultaneous design of proportional–integral–derivative controller and measurement filter by optimisation. Control Theory and Applications. 11(3), pp. 341-348. IET. 2017.
[4] K. Soltesz, A. Cervin. Is PID a good choice?. 3rd Conference on Advances in Proportional-Integral-Derivative Control. IFAC. Ghent, Belgium. 2018. Submitted manuscript.
[5] http://cvxr.com
[6] https://mosek.com
......@@ -44,6 +44,7 @@ opts = defaults(P,'optIAE',K2); % Hot start with filtered PID controller
Pc = P;
T = opts3.T; % Step response horizon
h = opts3.h; % Sampling period
t = 0:h:T; % Time vector to use in step response
wN = pi/h; % Nyquist frequency
w ={opts1.w(1),wN}; % Angular frequency san
......@@ -58,9 +59,17 @@ for design = 1:3
G = loopsens(P,K);
figure(1)
step(G.PSi,T)
y = step(G.PSi,t);
if design == 3
plot(t+h/2,y)
else
plot(t,y)
end
hold on
IAE = trapz(t,abs(y));
eval(['IAE' num2str(design) '=IAE'])
figure(2)
bode(G.So,w)
hold on
......@@ -76,12 +85,6 @@ for design = 1:3
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
......@@ -109,7 +112,4 @@ 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
title('Controller - K')
\ No newline at end of file
......@@ -162,4 +162,4 @@ 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
end
......@@ -61,7 +61,7 @@ K = K(p);
[e,t] = step(P*S(p));
[e,t] = resample(e,t);
J = trapz(t,abs(e));
% Gradient
e_p = step(-P^2*K_p*S(p)^2,t);
J_p = trapz(t,mdot(sign(e),e_p));
......@@ -79,7 +79,7 @@ K = K(p);
cS = [];
Sm_p = [];
end
% Complementary sensitivity
if Mt < inf
Tw = 1-Sw;
......@@ -90,7 +90,7 @@ K = K(p);
cT = [];
Tm_p = [];
end
% Noise sensitivity
if Mks < inf
Kw = fr(K(p));
......@@ -102,13 +102,13 @@ K = K(p);
cKS = [];
KSm_p = [];
end
% Constraints
c = [cS cT cKS];
cEq = [];
% Gradients
c_p = [Sm_p;Tm_p;KSm_p].';
cEq_p = [];
end
end
\ No newline at end of file
end
......@@ -65,7 +65,7 @@ for iter = 1:maxIter
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)-prevObj) < stopTol*prevObj
break;
......@@ -76,4 +76,4 @@ end
s = zpk('s');
K1 = [1 1/s s];
K1 = K1(1:n);
K = K1*p;
\ No newline at end of file
K = K1*p;
......@@ -106,7 +106,7 @@ K = minreal(K(p));
cS = [];
Sm_p = [];
end
% Complementary sensitivity
if Mt < inf
Tw = 1-Sw;
......@@ -129,13 +129,13 @@ K = minreal(K(p));
cKS = [];
KSm_p = [];
end
% Constraints
c = [cS cT cKS];
cEq = [];
% Gradients
c_p = [Sm_p;Tm_p;KSm_p].';
cEq_p = [];
end
end
\ No newline at end of file
end
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