Pole Placement Techniques Part 2

My task is that i should

1. Determine the system's transfer function and express the system's state equations in phase-variable form.
2. Identify a set of state-feedback gains using the pole placement technique to achieve specific performance criteria.
3. Convert the obtained set of state-feedback gains in the phase-variable system to the original system representation.
4. Verify that the set of gains in step (iii) successfully places the closed-loop poles at the desired positions.

This is my code thus far and my answer was deemed incompleted to which i've asked my peers and we all arrived at the same conclusion regarding task 3 and 4. Can you point me to the right direction on how i should approach? Thank you in advance.

clc; clear

% (1)
% Define the state-space matrices
A = [0 -83.33; 500 -10];    
B = [166.67; 0];
C = [0 1];
D = 0;

% Create a state-space system
sys = ss(A, B, C, D)

sys =
 
  A = 
           x1      x2
   x1       0  -83.33
   x2     500     -10
 
  B = 
          u1
   x1  166.7
   x2      0
 
  C = 
       x1  x2
   y1   0   1
 
  D = 
       u1
   y1   0
 
Continuous-time state-space model.

% Convert the state-space system to transfer function
disp('The system transfer function:');

The system transfer function:

tf_sys = tf(sys)

tf_sys =
 
          83335
  ---------------------
  s^2 + 10 s + 4.167e04
 
Continuous-time transfer function.

% Convert to phase-variable form
[Apv, Bpv, Cpv, Dpv] = canon(A, B, C, D, 'companion');

% Display the state-space matrices in phase-variable form
disp('Phase-Variable Form Matrices:');

Phase-Variable Form Matrices:

disp('Apv:');

Apv:

disp(Apv);

           0      -41665
           1         -10

disp('Bpv:');

Bpv:

disp(Bpv);

     1
     0

disp('Cpv:');

Cpv:

disp(Cpv);

           0       83335

disp('Dpv:');

Dpv:

disp(Dpv);

     0

%(2)
% Specify desired 20% overshoot and 0.5s settling time
OS = 0.2; %overshoot
ST = 0.5; %settling time

% Calculate natural frequency and damping ratio of system
zeta = -log(OS) / sqrt(pi^2 + (log(OS))^2);
omega_n = 4 / (zeta * ST);
disp('omega_n=');

omega_n=

disp(omega_n);

   17.5458

disp('zeta=');

zeta=

disp(zeta);

    0.4559

% Calculate desired poles
sigma = -zeta * omega_n;
omega_d = omega_n * sqrt(1 - zeta^2);
desiredPoles = [sigma + 1i*omega_d, sigma - 1i*omega_d];
disp('The desired poles for the system :')

The desired poles for the system :

disp(desiredPoles);

  -8.0000 +15.6159i  -8.0000 -15.6159i

% Use the place command to find state-feedback gains
K = place(Apv, Bpv, desiredPoles);
disp('State-feedback gains (K):');

State-feedback gains (K):

disp(K);

   1.0e+04 *

    0.0006   -4.1417

sys_cl = ss(Apv - Bpv*K, Bpv, Cpv, Dpv); % Create the closed-loop state-space model
N = dcgain(sys_cl) %Steady-state value

N = 270.6958

Sys_cl = ss(Apv - Bpv*K, Bpv/N, Cpv, Dpv);
%step(sys_cl); % Plot the step response of the closed-loop system 

%(3)
% Diagonalize matrix A
[T, ~] = eig(A);
%T = V; 
%D = 0;

% Display transformation matrix T
disp('Transformation matrix (T):');

Transformation matrix (T):

disp(T);

   0.0093 + 0.3778i   0.0093 - 0.3778i
   0.9258 + 0.0000i   0.9258 + 0.0000i

% Calculate state-feedback gains in the original system
K_original = K / T;
disp('State-feedback gains in the original system (K_original):');

State-feedback gains in the original system (K_original):

disp(K_original)

   1.0e+04 *

   0.0000 - 5.4815i  -2.2365 + 0.0548i

%(4)
% Calculate closed-loop system matrix A_cl
A_cl = Apv - Bpv * K_original;
disp('The closed-loop system matrix A_cl:')

The closed-loop system matrix A_cl:

disp(A_cl);

   1.0e+04 *

   0.0000 + 5.4815i  -1.9300 - 0.0548i
   0.0001 + 0.0000i  -0.0010 + 0.0000i

sys_cl2 = ss(A_cl, Bpv, Cpv, Dpv)

sys_cl2 =
 
  A = 
                   x1              x2
   x1     0+5.48e+04i  -1.93e+04-548i
   x2               1             -10
 
  B = 
       u1
   x1   1
   x2   0
 
  C = 
              x1         x2
   y1          0  8.334e+04
 
  D = 
       u1
   y1   0
 
Continuous-time state-space model.

% Calculate closed-loop poles
closedLoopPoles = eig(A_cl);

% Display the results
disp('Desired pole locations:');

Desired pole locations:

disp(desiredPoles);

  -8.0000 +15.6159i  -8.0000 -15.6159i

disp('Actual closed-loop pole locations:');

Actual closed-loop pole locations:

disp(closedLoopPoles);

   1.0e+04 *

  -0.0000 + 5.4815i
  -0.0010 - 0.0000i

% Check if the poles match the desired locations
if isequal(round(closedLoopPoles, 4), round(desiredPoles, 4))
    disp('The closed-loop poles match the desired locations.');
else
    disp('The closed-loop poles do not match the desired locations.');
end

The closed-loop poles do not match the desired locations.

%Plot Step Response
step(sys), hold on
step(Sys_cl), hold on
step(sys_cl2), grid on

Warning: The data cannot be plotted because it is not real valued.

legend('Original System', 'Compensated System', 'TBH IDK')

Answer 1

[Paul]

Hi MUHAMMAD,

See some comments interspersed with the code below.

clc; clear

Here, we define the sysem in original form.

% (1)
% Define the state-space matrices
A = [0 -83.33; 500 -10];    
B = [166.67; 0];
C = [0 1];
D = 0;
% Create a state-space system
sys = ss(A, B, C, D);

% Convert the state-space system to transfer function
tf_sys = tf(sys);

This form of the use of canon is obsolete. Consider changing to canon, but even that function is no longer recommended. I believe the current function is compreal. But there's nothing fundamentally wrong here.

As @Sam Chak pointed out in your other question, this form might not be phase-variable form, but we will assume it is for this discussion.

% Convert to phase-variable form
[Apv, Bpv, Cpv, Dpv] = canon(A, B, C, D, 'companion');


%(2)
% Specify desired 20% overshoot and 0.5s settling time
OS = 0.2; %overshoot
ST = 0.5; %settling time
% Calculate natural frequency and damping ratio of system
zeta = -log(OS) / sqrt(pi^2 + (log(OS))^2);
omega_n = 4 / (zeta * ST);

% Calculate desired poles
sigma = -zeta * omega_n;
omega_d = omega_n * sqrt(1 - zeta^2);
desiredPoles = [sigma + 1i*omega_d, sigma - 1i*omega_d];

% Use the place command to find state-feedback gains
K = place(Apv, Bpv, desiredPoles);

sys_cl = ss(Apv - Bpv*K, Bpv, Cpv, Dpv); % Create the closed-loop state-space model
N = dcgain(sys_cl); %Steady-state value

Sys_cl = ss(Apv - Bpv*K, Bpv/N, Cpv, Dpv);

Verify that the desired poles match the closed loop poles in phase-variable form.

desiredPoles

desiredPoles = 
  -8.0000 +15.6159i  -8.0000 -15.6159i

pole(Sys_cl).'

ans = 
  -8.0000 +15.6159i  -8.0000 -15.6159i

So far, so good, at least as far as placing the poles in the phase-variable form.

This next step isn't correct. The matrix T would be used to as the transformation between original form and modal form. But you need the transformation from original form to phase-variable form.

%(3)
% Diagonalize matrix A
[T, ~] = eig(A);

If T was the correct transformation, this next line might be correct. Depending on how T is defined, it might have to be changed as in the commented line

% Calculate state-feedback gains in the original system
K_original = K / T;
%K_original = K * T;

If K_original is correct, this next step is still incorrect. The feedback gain matrix should be applied to the system in its original form, not in phase-variable form.

%(4)
% Calculate closed-loop system matrix A_cl
A_cl = Apv - Bpv * K_original;

sys_cl2 = ss(A_cl, Bpv, Cpv, Dpv);

Summary:

Get the correct transformation T. Hint: check for an additional ouput from canon() (or compreal() if you decide to use that) if the form returned from canon is the form you want as phase-variable form. Otherwise, you'll have to properly define what phase-variable form means and then find its transformation T from original form.

Use the correct system in step (4).

[MUHAMMAD ANAS SULHI]

Thank you for the advice, Mr @Paul and @Sam Chak. After a lot of testing and a lot of reading, apparently the 'companion' can be used to find Transformation Matrix T, though i dont know if i implemented it correctly but i've get a plausible answer for Task 3 and Task 4 where the closed-loop poles match the desired locations.

% Define the state-space matrices
A = [0 -83.33; 500 -10];    
B = [166.67; 0];
C = [0 1];
D = 0;

% Create a state-space system
sys = ss(A, B, C, D)

sys =
 
  A = 
           x1      x2
   x1       0  -83.33
   x2     500     -10
 
  B = 
          u1
   x1  166.7
   x2      0
 
  C = 
       x1  x2
   y1   0   1
 
  D = 
       u1
   y1   0
 
Continuous-time state-space model.

% Convert the state-space system to transfer function
disp('The system transfer function:');

The system transfer function:

tf_sys = tf(sys)

tf_sys =
 
          83335
  ---------------------
  s^2 + 10 s + 4.167e04
 
Continuous-time transfer function.

% Convert to phase-variable form
[Apv, Bpv, Cpv, Dpv, T] = canon(A, B, C, D, 'companion');
spv = ss(Apv',Cpv',Bpv',Dpv')

spv =
 
  A = 
               x1          x2
   x1           0           1
   x2  -4.166e+04         -10
 
  B = 
              u1
   x1          0
   x2  8.334e+04
 
  C = 
       x1  x2
   y1   1   0
 
  D = 
       u1
   y1   0
 
Continuous-time state-space model.

% Specify desired 20% overshoot and 0.5s settling time
OS = 0.2; %overshoot
ST = 0.5; %settling time

% Calculate natural frequency and damping ratio of system
zeta = -log(OS) / sqrt(pi^2 + (log(OS))^2);
omega_n = 4 / (zeta * ST);
disp('omega_n=');

omega_n=

disp(omega_n);

   17.5458

disp('zeta=');

zeta=

disp(zeta);

    0.4559

% Calculate desired poles
sigma = -zeta * omega_n;
omega_d = omega_n * sqrt(1 - zeta^2);
desiredPoles = [sigma + 1i*omega_d, sigma - 1i*omega_d];
disp('The desired poles for the system :')

The desired poles for the system :

disp(desiredPoles);

  -8.0000 +15.6159i  -8.0000 -15.6159i

% Use the place command to find state-feedback gains
K = place(Apv, Bpv, desiredPoles);
disp('State-feedback gains (K):');

State-feedback gains (K):

disp(K);

   1.0e+04 *

    0.0006   -4.1417

% Create the closed-loop phase-variable system
sys_cl = ss(Apv - Bpv*K, Bpv, Cpv, Dpv);
N1 = dcgain(sys_cl) %Steady-state value

N1 = 270.6958

Sys_cl = ss(Apv - Bpv*K, Bpv/N1, Cpv, Dpv);

% Display transformation matrix T
disp('Transformation matrix (T):');

Transformation matrix (T):

T_inv = inv(T)

T_inv = 2×2

1.0e+04 * 0.0167 0 0 8.3335

% Calculate state-feedback gains in the original system
K_original = K / T_inv;
disp('State-feedback gains in the original system (K_original):');

State-feedback gains in the original system (K_original):

disp(K_original)

    0.0360   -0.4970

% Calculate closed-loop system matrix A_cl
A_cl = A - B * K_original;
disp('The closed-loop system matrix A_cl:')

The closed-loop system matrix A_cl:

disp(A_cl);

   -6.0000   -0.4957
  500.0000  -10.0000

% Create the closed-loop original system
sys_cl2 = ss(A_cl, B, C, D)

sys_cl2 =
 
  A = 
            x1       x2
   x1       -6  -0.4957
   x2      500      -10
 
  B = 
          u1
   x1  166.7
   x2      0
 
  C = 
       x1  x2
   y1   0   1
 
  D = 
       u1
   y1   0
 
Continuous-time state-space model.

N2 = dcgain(sys_cl2)

N2 = 270.6958

Sys_cl2 = ss(A_cl, B/N2, C, D)

Sys_cl2 =
 
  A = 
            x1       x2
   x1       -6  -0.4957
   x2      500      -10
 
  B = 
           u1
   x1  0.6157
   x2       0
 
  C = 
       x1  x2
   y1   0   1
 
  D = 
       u1
   y1   0
 
Continuous-time state-space model.

% Calculate closed-loop poles
closedLoopPoles = eig(A_cl);

% Display the results
disp('Desired pole locations:');

Desired pole locations:

disp(desiredPoles);

  -8.0000 +15.6159i  -8.0000 -15.6159i

disp('Actual closed-loop pole locations:');

Actual closed-loop pole locations:

disp(closedLoopPoles);

  -8.0000 +15.6159i
  -8.0000 -15.6159i

figure (1)
step(sys), hold on
step(Sys_cl), grid on
legend('State-Space System', 'Phase-Variable System')

figure (2)
step(sys), hold on
step(Sys_cl2), grid on
legend('State-Space System', 'Original System')

[Sam Chak]

Hi @MUHAMMAD ANAS SULHI,

You're catching on quickly. Even though the question might not explicitly state it, I'm confident that you can create a table comparing the step response characteristics of the Original system and the Compensated system. Feel free to include remarks in the table, highlighting the percentage of improvement. Given that these are linear systems, you can also utilize inverse Laplace transforms ilaplace() to derive the analytical equations for both systems.

In your final remarks of the assignment, share your insights on what you've learned from this exercise and delve into the significance of the transformation matrix in the control design of real systems. Keep up the great work! (By the way, the legends can be rephrased to convey clearer descriptions of the comparison).

% Define the state-space matrices
A = [0 -83.33; 500 -10];    
B = [166.67; 0];
C = [0 1];
D = 0;

% Create a state-space system
sys = ss(A, B, C, D);

% Convert to phase-variable form
[Apv, Bpv, Cpv, Dpv, T] = canon(A, B, C, D, 'companion');
spv = ss(Apv',Cpv',Bpv',Dpv');

% Specify desired 20% overshoot and 0.5s settling time
OS = 0.2; %overshoot
ST = 0.5; %settling time

% Calculate natural frequency and damping ratio of system
zeta = -log(OS) / sqrt(pi^2 + (log(OS))^2);
omega_n = 4 / (zeta * ST);

% Calculate desired poles
sigma = -zeta * omega_n;
omega_d = omega_n * sqrt(1 - zeta^2);
desiredPoles = [sigma + 1i*omega_d, sigma - 1i*omega_d];

% Use the place command to find state-feedback gains
K = place(Apv, Bpv, desiredPoles);

% Create the closed-loop phase-variable system
sys_cl = ss(Apv - Bpv*K, Bpv, Cpv, Dpv);
N1 = dcgain(sys_cl); % Steady-state value
Sys_cl = ss(Apv - Bpv*K, Bpv/N1, Cpv, Dpv);

% Display transformation matrix T
% disp('Transformation matrix (T):');
T_inv = inv(T);
% Calculate state-feedback gains in the original system
K_original = K / T_inv;
% disp('State-feedback gains in the original system (K_original):');
% disp(K_original)

% Calculate closed-loop system matrix A_cl
A_cl = A - B*K_original;

% Create the closed-loop original system
sys_cl2 = ss(A_cl, B, C, D);
N2 = dcgain(sys_cl2);
Sys_cl2 = ss(A_cl, B/N2, C, D);

stepinfo(sys)

ans = struct with fields:
         RiseTime: 0.0052
    TransientTime: 0.7714
     SettlingTime: 0.7714
      SettlingMin: 0.2854
      SettlingMax: 3.8521
        Overshoot: 92.5910
       Undershoot: 0
             Peak: 3.8521
         PeakTime: 0.0154

stepinfo(Sys_cl2)

ans = struct with fields:
         RiseTime: 0.0888
    TransientTime: 0.4746
     SettlingTime: 0.4746
      SettlingMin: 0.9334
      SettlingMax: 1.2000
        Overshoot: 19.9997
       Undershoot: 0
             Peak: 1.2000
         PeakTime: 0.2015

step(sys), hold on
step(Sys_cl2), grid on
legend('Original System', 'Compensated System')

[Paul]

@MUHAMMAD ANAS SULHI

Looks like you have it all sorted out, assuming that the companion form is acceptable as "phase-variable form." If it's not, you'll have to do some additional reading and testing to determine how to transform your system, i.e., find the T, from its origingal form to phase-variable form, which is also called controllable canonical form or controller canonical form.

These lines are basically undoing each other. The first line finds the inverse of T and the second effectively multiplies K by the inverse of the inverse of T, which is just T

T_inv = inv(T)
% Calculate state-feedback gains in the original system
K_original = K / T_inv;

So it would be more efficient and accurate to just do one line

K_original = K * T;

though it probably doesn't matter in this problem.

When you form the closed loop system, such as

sys_cl2 = ss(A_cl, B, C, D)

don't forget the possbility that D could be non-zero. It is zero in this problem, but if it weren't you'd have to make an adjustment to the line above to the the correct form of sys_cl2.C.

With the same idea, the input to the sysc_cl2 flows through the D matrix as well, so any dc gain adjustment, as being done here

Sys_cl2 = ss(A_cl, B/N2, C, D)

would also have to scale D.