# How to design a robust PI controller based on the Sensitivity peak value

2 ビュー (過去 30 日間)
Rinitha Rajan 2023 年 7 月 27 日
コメント済み: Rinitha Rajan 2023 年 7 月 31 日
I've been trying to implement this paper in MATLAB and here is the code so far.
% Parameters
f = 2000; %Hz
Tm = 1e-4;
Tr = 2.5e-5;
T1 = 0.0041;
T2 = 1.0024e-4;
Kr =1.56;
Kc = 16.3333;
Tc = 1.0024e-4;
K = T1/(2*Tr);
W = 2*pi*f;
Wn = sqrt((K+1)/(T1*Tr));
num_zeta = (T1+Tr)/(T1*Tr);
den_zeta = 2*(sqrt((K+1)/(T1*Tr)));
zeta = num_zeta/den_zeta;
% Gp
A = 9.6767e+08 -1i*5.0572e+08;
B = 1;
Bdot = 0;
D = A^2 + B^2;
%sensitivity analysis equation:
% solving for M
syms M;
M_eqn = D*alpha^2*(beta^2+W^(-2)) + 2*alpha*A*beta + 2*alpha*beta*(W^(-1)) + 1 - M^(-2) >= 0;
Error using alpha
Too many output arguments.
M_val = solve(M_eqn, M);
num_alpha = -2*Adot*beta - 2*Bdot*W^(-1) + 2*B*W^(-2);
den_alpha = 2*E*beta^(2) + 2*E*W^(-2) - 2*D*W^(-3);
alpha = num_alpha/den_alpha;
x2 = 4*B^(2)*D*W^(-4) + 8*A*Adot*D*W^(-3) - 8*B*Bdot*D*W^(-3) + 8*B^(2)*E*W^(-3) - 8*D*E*W^(-3) + 8*D*E*W^(-3)*M^(-2) + 4*Adot^(2)*D*W^(-2) + 4*Bdot^(2)*D*W^(-2) - 8*A*Adot*E*W^(-2) - 8*B*Bdot*E*W^(-2) + 8*E^(2)*W^(-2) - 8*E^(-2)*M^(-2);
x1 = -8*A*B*D*W^(-5) + 8*A*Bdot*D*W^(-4) + 8*A*B*E*W^(-4) + 8*Adot*Bdot*D*W^(-3) - 8*Adot*B*E*W^(-3) - 8*A*Bdot*E*W^(-3);
x0 = -4*B^(2)*D*W^(-6) + 4*D^(2)*W^(-6) - 4*D^(2)*W^(-6)*M^(-2) + 8*B^(2)*E*W^(-5) - 8*D*E*W^(-5) + 8*D*E*W^(-5)*M^(-2) * 4*B^(2)*D*W^(-4) - 8*B*Bdot*E*W^(-4) + 4*E^(2)*W^(-4) - 4*E^(2)*W^(-4)*M^(-2);
beta_val = vpasolve(x4*beta^(4)*x3*beta^(3) + x2*beta^(2) + x1*beta + x0 == 0, beta);
(The parameters have been given to me.)
I have been getting this error: "Error using alpha. Too many output arguments.
I am not as well-versed with MALTAB as I'd like to be..so please bear with me!
And if possible, could I also know if I'm on the right track and what I should be doing further from here?
##### 4 件のコメント2 件の古いコメントを表示2 件の古いコメントを非表示
Sam Chak 2023 年 7 月 28 日
I suggest to update / change the title of your Question to
"How to design a robust PI controller based on the Sensitivity peak value ?"
People who are good at the frequency response method and solving symbolic equations should be able to guide you.
Rinitha Rajan 2023 年 7 月 31 日
Thank you for the suggestion!

サインインしてコメントする。

### 回答 (3 件)

Steven Lord 2023 年 7 月 27 日
Nowhere have you defined a variable named alpha, so when MATLAB sees alpha in your code it calls the alpha function with one output argument. It will use that output argument in your equation. But the alpha function doesn't return any output arguments so the call to the function throws an error.
Define your alpha variable before you use it.
##### 1 件のコメント-1 件の古いコメントを表示-1 件の古いコメントを非表示
Rinitha Rajan 2023 年 7 月 27 日
Thank you so much! I did not know that alpha and beta are inbuilt functions..!

サインインしてコメントする。

VBBV 2023 年 7 月 27 日
Define or write these lines of code before the symbolic equation M_eqn. This equation contains two undefined variables, alpha and beta
For variable alpha, It seems you have written after the equation M_eqn,
% place these lines before
num_alpha = -2*Adot*beta - 2*Bdot*W^(-1) + 2*B*W^(-2);
den_alpha = 2*E*beta^(2) + 2*E*W^(-2) - 2*D*W^(-3);
% this variable alpha
alpha = num_alpha/den_alpha;
for variable beta give some value say, 2; and then solve the equation M_eqn
% Parameters
f = 2000; %Hz
Tm = 1e-4;
Tr = 2.5e-5;
T1 = 0.0041;
T2 = 1.0024e-4;
Kr =1.56;
Kc = 16.3333;
Tc = 1.0024e-4;
K = T1/(2*Tr);
W = 2*pi*f;
Wn = sqrt((K+1)/(T1*Tr));
num_zeta = (T1+Tr)/(T1*Tr);
den_zeta = 2*(sqrt((K+1)/(T1*Tr)));
zeta = num_zeta/den_zeta;
% Gp
A = 9.6767e+08 -1i*5.0572e+08;
B = 1;
Bdot = 0;
D = A^2 + B^2;
% define beta also
beta = 2;
% define them here
num_alpha = -2*Adot*beta - 2*Bdot*W^(-1) + 2*B*W^(-2);
den_alpha = 2*E*beta^(2) + 2*E*W^(-2) - 2*D*W^(-3);
alpha = num_alpha/den_alpha;
%sensitivity analysis equation:
% solving for M
syms M;
M_eqn = D*alpha^2*(beta^2+W^(-2)) + 2*alpha*A*beta + 2*alpha*beta*(W^(-1)) + 1 - M^(-2) >= 0;
M_val = solve(M_eqn, M);
x2 = 4*B^(2)*D*W^(-4) + 8*A*Adot*D*W^(-3) - 8*B*Bdot*D*W^(-3) + 8*B^(2)*E*W^(-3) - 8*D*E*W^(-3) + 8*D*E*W^(-3)*M^(-2) + 4*Adot^(2)*D*W^(-2) + 4*Bdot^(2)*D*W^(-2) - 8*A*Adot*E*W^(-2) - 8*B*Bdot*E*W^(-2) + 8*E^(2)*W^(-2) - 8*E^(-2)*M^(-2);
x1 = -8*A*B*D*W^(-5) + 8*A*Bdot*D*W^(-4) + 8*A*B*E*W^(-4) + 8*Adot*Bdot*D*W^(-3) - 8*Adot*B*E*W^(-3) - 8*A*Bdot*E*W^(-3);
x0 = -4*B^(2)*D*W^(-6) + 4*D^(2)*W^(-6) - 4*D^(2)*W^(-6)*M^(-2) + 8*B^(2)*E*W^(-5) - 8*D*E*W^(-5) + 8*D*E*W^(-5)*M^(-2) * 4*B^(2)*D*W^(-4) - 8*B*Bdot*E*W^(-4) + 4*E^(2)*W^(-4) - 4*E^(2)*W^(-4)*M^(-2);
beta_val = vpasolve(x4*beta^(4)*x3*beta^(3) + x2*beta^(2) + x1*beta + x0 == 0, beta)
beta_val =
##### 4 件のコメント2 件の古いコメントを表示2 件の古いコメントを非表示
VBBV 2023 年 7 月 27 日
Ok, Then define beta as symbolic variable
% Parameters
f = 2000; %Hz
Tm = 1e-4;
Tr = 2.5e-5;
T1 = 0.0041;
T2 = 1.0024e-4;
Kr =1.56;
Kc = 16.3333;
Tc = 1.0024e-4;
K = T1/(2*Tr);
W = 2*pi*f;
Wn = sqrt((K+1)/(T1*Tr));
num_zeta = (T1+Tr)/(T1*Tr);
den_zeta = 2*(sqrt((K+1)/(T1*Tr)));
zeta = num_zeta/den_zeta;
% Gp
A = 9.6767e+08 -1i*5.0572e+08;
B = 1;
Bdot = 0;
D = A^2 + B^2;
% define beta as symbolic variable
syms beta
% define them here
num_alpha = -2*Adot*beta - 2*Bdot*W^(-1) + 2*B*W^(-2);
den_alpha = 2*E*beta^(2) + 2*E*W^(-2) - 2*D*W^(-3);
alpha = num_alpha/den_alpha;
%sensitivity analysis equation:
% solving for M
syms M beta;
M_eqn = D*alpha^2*(beta^2+W^(-2)) + 2*alpha*A*beta + 2*alpha*beta*(W^(-1)) + 1 - M^(-2) >= 0;
M_val = vpasolve(M_eqn, M);
x2 = 4*B^(2)*D*W^(-4) + 8*A*Adot*D*W^(-3) - 8*B*Bdot*D*W^(-3) + 8*B^(2)*E*W^(-3) - 8*D*E*W^(-3) + 8*D*E*W^(-3)*M^(-2) + 4*Adot^(2)*D*W^(-2) + 4*Bdot^(2)*D*W^(-2) - 8*A*Adot*E*W^(-2) - 8*B*Bdot*E*W^(-2) + 8*E^(2)*W^(-2) - 8*E^(-2)*M^(-2);
x1 = -8*A*B*D*W^(-5) + 8*A*Bdot*D*W^(-4) + 8*A*B*E*W^(-4) + 8*Adot*Bdot*D*W^(-3) - 8*Adot*B*E*W^(-3) - 8*A*Bdot*E*W^(-3);
x0 = -4*B^(2)*D*W^(-6) + 4*D^(2)*W^(-6) - 4*D^(2)*W^(-6)*M^(-2) + 8*B^(2)*E*W^(-5) - 8*D*E*W^(-5) + 8*D*E*W^(-5)*M^(-2) * 4*B^(2)*D*W^(-4) - 8*B*Bdot*E*W^(-4) + 4*E^(2)*W^(-4) - 4*E^(2)*W^(-4)*M^(-2);
beta_val = vpasolve(x4*beta^(4)*x3*beta^(3) + x2*beta^(2) + x1*beta + x0 == 0, beta)
beta_val =
Rinitha Rajan 2023 年 7 月 28 日
Thank you so much! Could I know how I can get multiple values for M? As there is a parametric plot later on, which is Kp vs Ki for different values of M..
As with vpasolve, M-val gives me "Empty sym: 0-by-1"

サインインしてコメントする。

Sam Chak 2023 年 7 月 28 日
From the excerpt below, the α and β are the parameters to be determined to make up for the Proportional–Integral controller that is given by
.
Thus,
and ,
and I don't think they are some kind of inbuilt functions.
Excerpt 1:
I think computing α and β is fairly straightforward as outlined in the paper (see Excerpt 2 below) so long as you properly define the real values for . The parameter α is calculated directly from the formula Eq. (7), but the parameter β is obtained by solving the 4th-degree polynomial equation (quartic equation) as shown in Eq. (8).
Excerpt 2:
If you can follow the framework outlined by the paper that involves mostly algebraic substitutions on the denominator of the Sensitivity function in Eq. (5), then you should be able to find α and β that satisfy the desired sensitivity inequality.
##### 3 件のコメント1 件の古いコメントを表示1 件の古いコメントを非表示
Rinitha Rajan 2023 年 7 月 28 日

And yes this does work, I'm not sure how you found the M value without solving the sensitivity equation?
I also tried this with the new transfer function (as I have mentioned above)! Unfornately alpha, beta, Kp and Ki are turning out to be complex..
Sam Chak 2023 年 7 月 28 日
I didn't find M. In fact, I just randomly assigned a value for M, so that I can have a complete 4th-degree polynomial equation (see Eq. (8)) to be solved to get the real value of β.
Once β is determined, I use the formula in Eq. (7) to calculate the value of α.
According to the design procedure in Perng et al., 2016, it depends on the desired Sensitivity value you set for M. Thus, you don't need to symbolically solve for M.
Excerpt:
In the following example, the PI controller is directly tuned so that it gives a phase margin of 60°. Can you find the Sensitivity value M here?
% Parameters
Tm = 1e-4;
Tr = 2.5e-5;
T1 = 0.0041;
T2 = 1.0024e-4;
Kr = 1.56;
% Plant transfer function Gp and its step response
s = tf('s');
Gp = (Kr*(Tm*s + 1))/((Tr*s + 1)*(T1*s + 1)*(T2*s + 1));
Gp = minreal(Gp)
Gp = 1.518e07 s + 1.518e11 ------------------------------------------ s^3 + 5.022e04 s^2 + 4.112e08 s + 9.733e10 Continuous-time transfer function.
step(Gp)
% Tune PI controller for Gp
[Gc, info] = pidtune(Gp, 'PI')
Gc = 1 Kp + Ki * --- s with Kp = 0.816, Ki = 588 Continuous-time PI controller in parallel form.
info = struct with fields:
Stable: 1 CrossoverFrequency: 492.9601 PhaseMargin: 60.0000
Gcl = feedback(Gc*Gp, 1)
Gcl = 1.24e07 s^2 + 1.329e11 s + 8.929e13 --------------------------------------------------------- s^4 + 5.022e04 s^3 + 4.236e08 s^2 + 2.302e11 s + 8.929e13 Continuous-time transfer function.
step(Gcl)
% Find Sensitivity function
X = AnalysisPoint('X');
T = feedback(Gp*X*Gc, 1);
S = getSensitivity(T, 'X')
Generalized continuous-time state-space model with 1 outputs, 1 inputs, 4 states, and the following blocks: X: Analysis point, 1 channels, 1 occurrences. Type "ss(S)" to see the current value and "S.Blocks" to interact with the blocks.
Gs = tf(S) % Sensitivity transfer function
Gs = From input "X" to output "X": s^4 + 5.022e04 s^3 + 4.112e08 s^2 + 9.733e10 s --------------------------------------------------------- s^4 + 5.022e04 s^3 + 4.236e08 s^2 + 2.302e11 s + 8.929e13 Continuous-time transfer function.
bode(Gs)

サインインしてコメントする。

### カテゴリ

Help Center および File ExchangeTime-Domain Analysis についてさらに検索

R2021b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by