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

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;
Adot = -2.5133e+04 + 4.024e+04i;
Bdot = 0;
D = A^2 + B^2;
E = A*(Adot) + B*(Bdot);

%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;

x4 = 4*Adot*D - 8*A*Adot*E + 4*E^2 - 4*E^(2)*M^(-2);
x3 = -8*Adot*B*D*W^(-2) +8*A*B*E*W^(-2) + 8*Adot*Bdot*D*W^(-1) - 8*Adot*B*E*W^(-1) - 8*A*Bdot*E*W^(-1);
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?

[Sam Chak]

Hi @Rinitha Rajan

Need some clarifications.

Item 1. In Eq. (4), the left-hand side implies that the plant is modeled as function of s. But the right-hand side is a complex number in the form , where A and B are real numbers. Can you show the transfer function of ?

Item 2. Also, does the linear system exist without its complex conjugate ?

Item 3. In your code, A is defined as a complex number, not a real number, as implied in Perng et al., 2016. How exactly do you compute , the time derivative of ?

Clarifying these are important because we want to know whether you are computing correctly or not. If it is wrong right from the beginning, then it casts doubt over your numerical solution of α and β.

Excerpt from: Perng, JW., Hsieh, SC., Ma, LS. et al. Design of robust PI control systems based on sensitivity analysis and genetic algorithms. Neural Comput & Applic 29, 913–923 (2018). https://doi.org/10.1007/s00521-016-2506-2

% Gp
A    =  9.6767e+08 - 5.0572e+08i;
B    =  1;
Adot = -2.5133e+04 + 4.024e+04i;
Bdot =  0;
D    =  A^2 + B^2

D = 6.8063e+17 - 9.7874e+17i

E    =  A*Adot + B*Bdot

E = -3.9703e+12 + 5.1649e+13i

[Rinitha Rajan]

Hi! Yes of course, I understand that I have given incomplete information.

1. This was my initial Gp(s)

(I apologize for the handwritten pictures; I am currently on my phone!)

And the above was what I used to write the given code.

2. For the above Gp(s), I'd assume the complex conjugate does exist.

3. Adot in this case, is actually differentiation of A with respect to omega(W). This, I have confirmed with one of the reference papers mentioned.

As of now, I have received a new plant or a new transfer function function:

And I've followed the same procedure as such. Althought the issuse of whether my calculations for alpha and beta still persists. I followed through with the Tan Method part of the paper (Section 2.2) as I already had the required values for computations.

Here are my calculations:

And here is the revised code for the above for reference:

close all;

% Parameters
f = 20000; %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;
s = 1i*W;


% Gp

A = (1+(1i*W*Tr))*(1+(1i*W*T1))*(1+(1i*W*T2));
B = Kr*(1+(1i*W*Tm));
Adot = 1i*(T1+T2+Tr)-2*1i*W*((T1*T2)+(T1*Tr)+(T2*Tr))-3*1i*W^(2)*(T1*T2*Tr);
Bdot = 1i*Kr*Tm;
D = A^2 + B^2;
E = A*(Adot) + B*(Bdot);

% define beta 
beta = 1;
 
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);


x4 = 4*Adot^(2)*D - 8*A*Adot*E + 4*E^2 - 4*E^(2)*M^(-2);
x3 = -8*Adot*B*D*W^(-2) +8*A*B*E*W^(-2) + 8*Adot*Bdot*D*W^(-1) - 8*Adot*B*E*W^(-1) - 8*A*Bdot*E*W^(-1);
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*Bdot^(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);


% Tan Method: o - odd part and e - even part

Be = Kr;
Bo = 1i*W*Kr*Tm;
Ae = 1+(1i*W)^(2)*((T1*T2)+(Tr*T2)+(T1*Tr));
Ao =  (1i*W)*(T1+T2+Tr)+(1i*W)^(3)*(T1*T2*Tr);

P1 = -W^(2)*Bo*(-W^(2));
P2 = Be*(-W^(2));
P3 = W*Be*(-W^(2));
P4 = W*Bo*(-W^(2));
P5 = W^(2)*Ao*(-W^(2)); 
P6 = -W*Ae*(-W^(2));


num_Gp = Be*(-W^(2))+1i*W*Bo*(-W^(2));
den_Gp = Ae*(-W^(2))+1i*W*Ao*(-W^(2));
Gp = tf(num_Gp,den_Gp);
Kp = (P5*P4-P6*P2)/(P1*P4-P2*P3);
Ki = (P6*P1-P5*P3)/(P1*P4-P2*P3);

Again, thank you for helping me and please do let me know if anything I wrote is unclear!

[Sam Chak]

Hi @Rinitha Rajan,

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]

@Sam Chak

Thank you for the suggestion!

Answer 1

[Steven Lord]

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.

[Rinitha Rajan]

Thank you so much! I did not know that alpha and beta are inbuilt functions..!

Answer 2

[VBBV]

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;
Adot = -2.5133e+04 + 4.024e+04i;
Bdot = 0;
D = A^2 + B^2;
E = A*(Adot) + B*(Bdot);

% 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);


x4 = 4*Adot*D - 8*A*Adot*E + 4*E^2 - 4*E^(2)*M^(-2);
x3 = -8*Adot*B*D*W^(-2) +8*A*B*E*W^(-2) + 8*Adot*Bdot*D*W^(-1) - 8*Adot*B*E*W^(-1) - 8*A*Bdot*E*W^(-1);
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 = 

[VBBV]

Note also that alpha & beta are standard Matlab functions , which when used without input arguments will return such errors

[Rinitha Rajan]

Thank you so much! However the calculation of alpha and beta plus all these equations are dependant on each other, so I'm not sure if the values I've got are correct or not....

[VBBV]

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;
Adot = -2.5133e+04 + 4.024e+04i;
Bdot = 0;
D = A^2 + B^2;
E = A*(Adot) + B*(Bdot);

% 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);


x4 = 4*Adot*D - 8*A*Adot*E + 4*E^2 - 4*E^(2)*M^(-2);
x3 = -8*Adot*B*D*W^(-2) +8*A*B*E*W^(-2) + 8*Adot*Bdot*D*W^(-1) - 8*Adot*B*E*W^(-1) - 8*A*Bdot*E*W^(-1);
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]

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"

Answer 3

[Sam Chak]

Hi @Rinitha Rajan

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.

[Sam Chak]

Can you test if this approach works for your problem?

% Parameters
f    = 2000;            % Hz
W    = 2*pi*f;

% Gp
A    =  9.6767e+08;     % real number
B    =  1;              % real number
Adot = -2.5133e+04;     % real number
Bdot =  0;              % real number
D    =  A^2 + B^2;
E    =  A*(Adot) + B*(Bdot); 
M    =  0.01;           % define the desired Sensitivity value

% Define beta as a symbolic variable 
syms beta

% Find beta by solving Eq. (8):
x4       =   4*Adot*D - 8*A*Adot*E + 4*E^2 - 4*E^(2)*M^(-2);
x3       = - 8*Adot*B*D*W^(-2) +8*A*B*E*W^(-2) + 8*Adot*Bdot*D*W^(-1) - 8*Adot*B*E*W^(-1) - 8*A*Bdot*E*W^(-1);
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 = 

% Find alpha using the fomula in Eq. (7):
beta      =  beta_val(1);     % choosing the real solution
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

alpha = 

% PI control gains
Ki = alpha

Ki = 

Kp = alpha*beta

Kp = 

[Rinitha Rajan]

Hi! I have replied to your previose reply!

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]

Hi @Rinitha Rajan

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)