LTI Systems: Qualitative Analysis of the Step Response
Table of Contents
Introduction
With a Little Help from the 'help' Command!
Step Response of First Order Systems
Strictly Proper Systems
The Step Response Characteristic Parameters
Not Strictly Proper Systems
The Step Response Characteristic Parameters
Step Response of Second Order Systems
Strictly Proper Systems - Complex Poles, no Zeros
The Step Response Characteristic Parameters
Not Strictly Proper Systems - Complex Poles, one Zeros
The Step Response Characteristic Parameters
Systems of Order greater than 2 - Dominant Poles Approximation
Example - A Real Dominant Pole
The Dominant Pole Approximation
Example 2 - A Real Dominant Pole and a Zero close to the Origin
The Dominant Pole Approximation
Example 3 - Complex Dominant Poles
The Dominant Poles Approximation
Final Example - A Set of Dominant Poles and Zeros
The Dominant Poles Approximation
Introduction
In this live script, we will illustrate with examples the qualitative behaviour of the step response for proper, asymptotically stable LTI systems of first, second and higher order, without or with zeros, varying the position of the poles or the position of the zeros if present. Furthermore, we will illustrate how to evaluate the characteristic parameters of the step response (i.e. rise time, settling time, etc.) using the functions available in MATLAB's Control System Toolbox.
clear
clc
With a Little Help from the 'help' Command!
Here is the documentation, available with the help command, related to the Control System Toolbox function we will use in this live script
help tf
help step
help RespConfig
help stepplot
help stepinfo
Step Response of First Order Systems
Strictly Proper Systems
Let's consider
Remeber:
Let's explore what happens when the pole 
moves, starting close to the origin in 
and moving away from the origin towards 
on the real axis
moves, starting close to the origin in and moving away from the origin towards on the real axisLet's choose some pole positions in the Left Half Plane of the complex frequency s (LHP for short).
Np = 7; % 7 different pole positions -- we have 7 available colors for plots
tau_p = logspace(-1,1, Np); % Np time constants tau, logarithmically spaced,
% starting from -1 e-1 and ending with -1 e+1
muG = 2+8*rand(1); % a positive value in the range [2, 10]
Pay attention on this function:
help stepinfo
We will use it to have an accurate estimation of the characteristic parameters of the step response and we'll compare the results with the approximate formulas in Eq. (+).
G1sN = struct('tf',cell(1), ...
'y_step', cell(1), ...
't_step', cell(1)); % we'll store each different transfer function
G1sStepINFO = cell(Np); % and the step response parameters
figure('Units','pixels', 'Position',[1,1,1280, 960]);
% let's create a blank figure, as wide as possible
currAx = gca; % retrieve the current axes handler
hold on;
Tfinal = 1.5*5*max(tau_p); % highest than the longest settling time
for n=Np:-1:1
% for each pole position
tau = tau_p(n);
currGs = tf(muG, [tau +1] ); % 1st order, NO zero
[Ystep, Tstep] = step(currGs, Tfinal); % plot the current step response
plot(currAx,Tstep, Ystep, 'LineWidth', 1.5);
G1sN(n).tf = currGs;
G1sN(n).y_step = Ystep;
G1sN(n).t_step = Tstep;
G1sStepINFO{n} = stepinfo(currGs, "SettlingTimeThreshold", 0.01);
end % for n
formatSpec = '%.2f'; % see 'doc num2str'
legend(['\tau = ', num2str(tau_p(7), formatSpec) ], ...
['\tau = ', num2str(tau_p(6), formatSpec) ], ...
['\tau = ', num2str(tau_p(5), formatSpec) ], ...
['\tau = ', num2str(tau_p(4), formatSpec) ], ...
['\tau = ', num2str(tau_p(3), formatSpec) ], ...
['\tau = ', num2str(tau_p(2), formatSpec) ], ...
['\tau = ', num2str(tau_p(1), formatSpec) ], ...
'Location', 'best', 'Fontsize', 14)
Notice: the system with 
has the slowest step response, whereas the system with 
is the fastest one.
has the slowest step response, whereas the system with is the fastest one. The Step Response Characteristic Parameters
Let's analyse the step response parameters, comparing the expressions in Eq. (+) with the estimations obtained by applying the function stepinfo( ).
What does the function stepinfo( ) provide?
disp(G1sStepINFO{1})
settlingTformula = zeros(Np,1); % initialize the array used to store
% the settling time values obtained
% using the formula
settlingT_stepinfo = zeros(size(settlingTformula));
for n = 1:Np
settlingTformula(n) = 4.6 * tau_p(n);
settlingT_stepinfo(n) = G1sStepINFO{n}.SettlingTime;
end % for n
settlingT_table = table(settlingTformula, settlingT_stepinfo, tau_p', ...
'VariableNames',...
{'formula 4.6 tau', 'estimated by stepinfo', 'tau values'});
disp(settlingT_table)
As you can note, the values are very close similar for the settling time of each step response.
Pay attention: this is the only scenario when both approaches give the same results. Any guess why?
Not Strictly Proper Systems
Let's consider
Let's choose a fixed pole position in the Left Half Plane of the complex frequency s (LHP for short), a static gain and explore what happens if the zero moves from LHP to RHP (Right Half Plane).
Np = 7; % 7 different pole positions -- we have 7 available colors for plots
T_z = [-fliplr(logspace(-1,1, (floor(Np/2)))), ...
0, ...
logspace(-1,1, (floor(Np/2)))]; % Np time constants T, logarithmically spaced,
% starting from -1e-1 and ending with -1e+1
tau_p = 1.5;
muG = 2+8*rand(1); % a positive value in the range [2, 10]
G1sZN = struct('tf',cell(1), ...
'y_step', cell(1), ...
't_step', cell(1)); % we'll store each different transfer function
G1sZStepINFO = cell(Np); % and the step response parameters
figure('Units','pixels', 'Position',[1,1,1280, 960]);
% let's create a blank figure, as wide as possible
currAx = gca; % retrieve the current axes handler
hold on;
Tfinal = 1.5*5*max(tau_p); % highest than the longest settling time
for n=1:Np
% for each pole position
Tz = T_z(n);
currGs = muG* tf([Tz 1], [tau_p +1] ); % 1st order, NO zero
[Ystep, Tstep] = step(currGs, Tfinal); % plot the current step response
plot(currAx,Tstep, Ystep, 'LineWidth', 1.5);
G1sZN(n).tf = currGs;
G1sZN(n).y_step = Ystep;
G1sZN(n).t_step = Tstep;
G1sZStepINFO{n} = stepinfo(currGs, "SettlingTimeThreshold", 0.01);
end % for n
formatSpec = '%.2f'; % see 'doc num2str'
legend(['T = ', num2str(T_z(1), formatSpec) ], ...
['T = ', num2str(T_z(2), formatSpec) ], ...
['T = ', num2str(T_z(3), formatSpec) ], ...
['T = ', num2str(T_z(4), formatSpec) ], ...
['T = ', num2str(T_z(5), formatSpec) ], ...
['T = ', num2str(T_z(6), formatSpec) ], ...
['T = ', num2str(T_z(7), formatSpec) ], ...
'Location', 'best', 'Fontsize', 14)
The Step Response Characteristic Parameters
Let's analyse the step response parameters, comparing the expressions in Eq. (+) with the estimations obtained by applying the function stepinfo( ).
settlingTformula = zeros(Np,1); % initialize the array used to store
% the settling time values obtained
% using the formula
settlingT_stepinfo = zeros(size(settlingTformula));
polePOS = -(1/tau_p)*ones(Np,1);
zeroPOS = zeros(Np,1);
for n = 1:Np
settlingTformula(n) = 4.6 * tau_p;
settlingT_stepinfo(n) = G1sZStepINFO{n}.SettlingTime;
zeroPOS(n) = 1/T_z(n);
end % for n
settlingT_table = table(settlingTformula, settlingT_stepinfo, ...
polePOS, zeroPOS, ...
'VariableNames',...
{'formula 4.6 tau', 'estimated by stepinfo', ...
'pole p', 'zero z'});
disp(settlingT_table)
Step Response of Second Order Systems
Strictly Proper Systems - Complex Poles, no Zeros
Let's consider
Let's choose a constant value for 
and μ, and explore what happens if the damping factor ξ moves from 0 to 1. Let's analyse the step response parameters, comparing the expressions in Eq. ( ⋆ ) with the estimations obtained by applying the function stepinfo( ).
and μ, and explore what happens if the damping factor ξ moves from 0 to 1. Let's analyse the step response parameters, comparing the expressions in Eq. ( ⋆ ) with the estimations obtained by applying the function stepinfo( ). Np = 7; % 7 different pole positions -- we have 7 available colors for plots
xi_val = linspace(1e-1,1-1e-1,7); % Np damping factor values
omega_n_p = 4.5;
muG = 2+8*rand(1); % a positive value in the range [2, 10]
G2sN = struct('tf',cell(1), ...
'y_step', cell(1), ...
't_step', cell(1)); % we'll store each different transfer function
G2sStepINFO = cell(Np,1); % and the step response parameters
Tfinal = 1.5*5*(1/omega_n_p/1e-1); % highest than the longest settling time
figure('Units','pixels', 'Position',[1,1,1280, 960]);
% let's create a blank figure, as wide as possible
currAx = gca; % retrieve the current axes handler
hold on;
for n=1:Np
% for each pole position
currXi = xi_val(n);
currGs = tf(muG, [(1/omega_n_p.^2) (2*currXi/omega_n_p) +1] ); % 2nd order, NO zero
[Ystep, Tstep] = step(currGs, Tfinal); % plot the current step response
plot(currAx,Tstep, Ystep, 'LineWidth', 1.5);
G2sN(n).tf = currGs;
G2sN(n).y_step = Ystep;
G2sN(n).t_step = Tstep;
G2sStepINFO{n} = stepinfo(currGs, "SettlingTimeThreshold", 0.01);
end % for n
legend(['\xi= ',num2str(xi_val(1), formatSpec) ], ...
['\xi= ', num2str(xi_val(2), formatSpec) ], ...
['\xi= ',num2str(xi_val(3), formatSpec) ], ...
['\xi= ',num2str(xi_val(4), formatSpec) ], ...
['\xi= ',num2str(xi_val(5), formatSpec) ], ...
['\xi= ',num2str(xi_val(6), formatSpec) ], ...
['\xi= ',num2str(xi_val(7), formatSpec) ],...
'Location', 'best', 'Fontsize', 14)
The Step Response Characteristic Parameters
Let's analyse the step response parameters, comparing the expressions in Eq. (+) with the estimations obtained by applying the function stepinfo( ).
settlingTformula = zeros(Np,1); % initialize the array used to store
% the settling time values obtained
% using the formula
settlingT_stepinfo = zeros(size(settlingTformula));
for n = 1:Np
settlingTformula(n) = 4.6 /(omega_n_p*xi_val(n));
settlingT_stepinfo(n) = G2sStepINFO{n}.SettlingTime;
end % for n
settlingT_table = table(settlingTformula, settlingT_stepinfo, xi_val', ...
'VariableNames',...
{'formula 4.6/(omega_n * xi)', 'estimated by stepinfo', 'damping factor'});
disp(settlingT_table)
Not Strictly Proper Systems - Complex Poles, one Zeros
Let's consider
Let's choose a constant value for , T and μ, and explore what happens if the damping factor ξ moves from 0 to 1. Let's analyse the step response parameters, comparing the expressions in Eq. ( ⋆ ) with the estimations obtained by applying the function stepinfo( ).
, T and μ, and explore what happens if the damping factor ξ moves from 0 to 1. Let's analyse the step response parameters, comparing the expressions in Eq. ( ⋆ ) with the estimations obtained by applying the function stepinfo( ). Np = 7; % 7 different pole positions -- we have 7 available colors for plots
xi_val = linspace(1e-1,1-1e-1,7); % Np damping factor values
omega_n_p = 4.5;
muG = 2+8*rand(1); % a positive value in the range [2, 10]
T = 1.25;
G2sN = struct('tf',cell(1), ...
'y_step', cell(1), ...
't_step', cell(1)); % we'll store each different transfer function
G2sStepINFO = cell(Np,1); % and the step response parameters
Tfinal = 1.5*5*(1/omega_n_p/1e-1); % highest than the longest settling time
figure('Units','pixels', 'Position',[1,1,1280, 960]);
% let's create a blank figure, as wide as possible
currAx = gca; % retrieve the current axes handler
hold on;
for n=1:Np
% for each pole position
currXi = xi_val(n);
currGs = muG*tf([T 1], [(1/omega_n_p.^2) (2*currXi/omega_n_p) +1] ); % 2nd order, NO zero
[Ystep, Tstep] = step(currGs, Tfinal); % plot the current step response
plot(currAx,Tstep, Ystep, 'LineWidth', 1.5);
G2sN(n).tf = currGs;
G2sN(n).y_step = Ystep;
G2sN(n).t_step = Tstep;
G2sStepINFO{n} = stepinfo(currGs, "SettlingTimeThreshold", 0.01);
end % for n
formatSpec = '%.2f'; % see 'doc num2str'
legend(['\xi = ',num2str(xi_val(1), formatSpec) ], ...
['\xi = ',num2str(xi_val(2), formatSpec) ], ...
['\xi = ',num2str(xi_val(3), formatSpec) ], ...
['\xi = ',num2str(xi_val(4), formatSpec) ], ...
['\xi = ',num2str(xi_val(5), formatSpec) ], ...
['\xi = ',num2str(xi_val(6), formatSpec) ], ...
['\xi = ',num2str(xi_val(7), formatSpec) ],...
'Location', 'best', 'Fontsize', 14)
The Step Response Characteristic Parameters
Let's analyse the step response parameters, comparing the expressions in Eq. (+) with the estimations obtained by applying the function stepinfo( ).
settlingTformula = zeros(Np,1); % initialize the array used to store
% the settling time values obtained
% using the formula
settlingT_stepinfo = zeros(size(settlingTformula));
for n = 1:Np
settlingTformula(n) = 4.6 /(omega_n_p*xi_val(n));
settlingT_stepinfo(n) = G2sStepINFO{n}.SettlingTime;
end % for n
settlingT_table = table(settlingTformula, settlingT_stepinfo, xi_val', ...
'VariableNames',...
{'formula 4.6/(omega_n * xi)', 'estimated by stepinfo', 'damping factor'});
disp(settlingT_table)
Systems of Order greater than 2 - Dominant Poles Approximation
When using the dominant poles approximation:
- It is essential to “preserve” the gain
- Zeros located close to the imaginary axis have to be adequately taken into account
This approximation is valid in qualitative analysis and for the initial and rough controller’s design steps. Refer to Part 7 of the course material.
Example - A Real Dominant Pole
Consider the LTI system described by the transfer function
clear
s = tf('s'); % the transfer function building helper
G1s = 250*(1+s/10)/( (1+20*s)*(1+s/5)*(1+(9/20)*s+s^2/16) ); % the transfer function
set(G1s,"Name", "G1");
G1_staticGAIN = dcgain(G1s) % the static gain of the LTI system
p_G1s = sort(pole(G1s), 'ascend') % the poles of the transfer function
zero(G1s) % the zeros of the transfer function
The poles and zeros configuration of the system is as follows
hf = figure('Units','pixels', 'Position',[1,1,1280, 960]);
hPZP = pzplot(hf, G1s);
grid on
get(hPZP)
get(hPZP.Responses(1))
% PZPlot properties:
% https://it.mathworks.com/help/ident/ref/controllib.chart.pzplot-properties.html
% https://it.mathworks.com/matlabcentral/answers/2173352-change-marker-colour-of-pzplot
%
hPZP.Responses(1).MarkerSize = 12;
hPZP.Responses(1).LineWidth = 2;
hPZP.Responses(1).Color = 'r';
get(hPZP.AxesStyle(1))
% see 'doc pzplot' & 'doc pzoptions'
hPZP.AxesStyle(1).FontSize = 14;
hPZP.AxesStyle(1).GridColor = [0,0,0];
hPZP.XLabel(1)
hPZP.XLabel(1).FontSize = 14;
hPZP.YLabel(1).FontSize = 14;
The pole closest to the origin is the dominant one.
The Dominant Pole Approximation
We may approximate the system 
using a very simple model: the static gain and the real dominant pole:
using a very simple model: the static gain and the real dominant pole:
tau1 = abs(1/p_G1s(1))
G1sApprox = tf([G1_staticGAIN], [tau1 1], "Name", "G_approx")
Let's compare the step response of system with the approximate model:
with the approximate model:hf = figure('Units','pixels', 'Position',[1,1,1280, 960]);
% let's create a blank figure
rpSP = stepplot(G1s, G1sApprox, [0, 180]);
rpSP
% refer to https://it.mathworks.com/help/control/ug/customizing-response-plots-from-the-command-line.html
rpSP.Characteristics(1).SettlingTime(1).Threshold = 0.01; % 1 %
rpSP.Characteristics(1).SettlingTime(1).Visible = "on";
rpSP.Characteristics(1).SteadyState(1).Visible = "on";
rpSP.LegendVisible = "on";
rpSP.LegendLocation = "best";
grid on;
rpSP.Responses
rpSP.Responses(1).LineWidth = 2;
rpSP.Responses(2).LineWidth = 2;
ylim([0, 260]);
rpSP.AxesStyle(1).FontSize=16;
rpSP.XLabel(1).FontSize = 16;
rpSP.YLabel(1).FontSize = 16;
G1info = stepinfo(G1s, "SettlingTimeThreshold", 0.01)
G1APPROXinfo = stepinfo(G1sApprox, "SettlingTimeThreshold", 0.01)
Example 2 - A Real Dominant Pole and a Zero close to the Origin
Consider the LTI system described by the transfer function
clear
s = tf('s'); % the transfer function building helper
G2s = 250*(1+50*s)/( (1+20*s)*(1+s/5)*(1+(9/20)*s+s^2/16) ); % the transfer function
set(G2s,"Name", "G2");
G2_staticGAIN = dcgain(G2s) % the static gain of the LTI system
p_G2s = sort(pole(G2s), 'ascend') % the poles of the transfer function
z_G2s = zero(G2s) % the zeros of the transfer function
The poles and zeros configuration of the system is as follows
hf = figure('Units','pixels', 'Position',[1,1,1280, 960]);
hPZP = pzplot(hf, G2s);
grid on
% get(hPZP)
% PZPlot properties:
% https://it.mathworks.com/help/ident/ref/controllib.chart.pzplot-properties.html
% https://it.mathworks.com/matlabcentral/answers/2173352-change-marker-colour-of-pzplot
%
hPZP.Responses(1).MarkerSize = 12;
hPZP.Responses(1).LineWidth = 2;
hPZP.Responses(1).Color = 'r';
% get(hPZP.AxesStyle(1))
% see 'doc pzplot' & 'doc pzoptions'
hPZP.AxesStyle(1).FontSize = 14;
hPZP.AxesStyle(1).GridColor = [0,0,0];
% hPZP.XLabel(1)
hPZP.XLabel(1).FontSize = 14;
hPZP.YLabel(1).FontSize = 14;
The pole closest to the origin is the dominant one. The zero close to the origin is also crucial to obtain a good apprimating model for the system.
The Dominant Pole Approximation
We may approximate the system 
using a very simple model: the static gain, the real dominant pole and the zero:
using a very simple model: the static gain, the real dominant pole and the zero:
tau2 = abs(1/p_G2s(1))
T2 = abs(1/z_G2s(1))
G2sApproxV1 = tf(G2_staticGAIN, [tau2 1], "Name", "G2_apprV1") % without the zero
G2sAPPROX = G2_staticGAIN * tf([T2 +1], [tau2 +1]) % with the zero
set(G2sAPPROX , "Name", "G2_apprV2")
Let's compare the step response of system with the approximate models:
with the approximate models:hf = figure('Units','pixels', 'Position',[1,1,1280, 960]);
% let's create a blank figure
rpSP2 =stepplot(G2s,G2sApproxV1, G2sAPPROX, [0, 180]);
grid on;
rpSP2
% refer to https://it.mathworks.com/help/control/ug/customizing-response-plots-from-the-command-line.html
rpSP2.Characteristics(1).SettlingTime(1).Threshold = 0.01; % 1 %
rpSP2.Characteristics(1).SettlingTime(1).Visible = "on";
rpSP2.Characteristics(1).SteadyState(1).Visible = "on";
rpSP2.LegendVisible = "on";
rpSP2.LegendLocation = "best";
grid on;
rpSP2.Responses
rpSP2.Responses(1).LineWidth = 2;
rpSP2.Responses(2).LineWidth = 2;
rpSP2.Responses(3).LineWidth = 2;
rpSP2.AxesStyle(1).FontSize=16;
rpSP2.XLabel(1).FontSize = 16;
rpSP2.YLabel(1).FontSize = 16;
G2info = stepinfo(G2s, "SettlingTimeThreshold", 0.01)
G2approx_V1_info = stepinfo(G2sApproxV1, "SettlingTimeThreshold", 0.01)
G2approx_info = stepinfo(G2sAPPROX, "SettlingTimeThreshold", 0.01)
Example 3 - Complex Dominant Poles
Consider the LTI system described by the transfer function
clear
s = tf('s'); % the transfer function building helper
G3s = 250*(1+(s/40))/( (1+(s/20))*(1+s/15)*(1+(9/20)*s+s^2/16) ); % the transfer function
set(G3s, "Name", "G3");
G3_staticGAIN = dcgain(G3s) % the static gain of the LTI system
p_G3s = sort(pole(G3s), 'ascend') % the poles of the transfer function
z_G3s = zero(G3s) % the zeros of the transfer function
The poles and zeros configuration of the system is as follows
hf = figure('Units','pixels', 'Position',[1,1,1280, 960]);
hPZP3 = pzplot(hf, G3s);
grid on
% get(hPZP)
% PZPlot properties:
% https://it.mathworks.com/help/ident/ref/controllib.chart.pzplot-properties.html
% https://it.mathworks.com/matlabcentral/answers/2173352-change-marker-colour-of-pzplot
%
hPZP3.Responses(1).MarkerSize = 12;
hPZP3.Responses(1).LineWidth = 2;
hPZP3.Responses(1).Color = 'r';
% get(hPZP.AxesStyle(1))
% see 'doc pzplot' & 'doc pzoptions'
hPZP3.AxesStyle(1).FontSize = 14;
hPZP3.AxesStyle(1).GridColor = [0,0,0];
% hPZP.XLabel(1)
hPZP3.XLabel(1).FontSize = 14;
hPZP3.YLabel(1).FontSize = 14;
The poles closest to the origin are the dominant ones.
The Dominant Poles Approximation
We may approximate the system 
using a very simple model: the static gain and the dominant poles
using a very simple model: the static gain and the dominant poles
G3sAPPROX = 250*(1)/( (1+(9/20)*s+s^2/16) ); % the approximating transfer function
set(G3sAPPROX, "Name", "G3_approx");
Let's compare the step response of system with the approximate model:
with the approximate model:hf = figure('Units','pixels', 'Position',[1,1,1280, 960]);
% let's create a blank figure
rpSP3 = stepplot(G3s, G3sAPPROX, [0, 180]);
rpSP3
% refer to https://it.mathworks.com/help/control/ug/customizing-response-plots-from-the-command-line.html
rpSP3.Characteristics(1).SettlingTime(1).Threshold = 0.01; % 1 %
rpSP3.Characteristics(1).SettlingTime(1).Visible = "on";
rpSP3.Characteristics(1).SteadyState(1).Visible = "on";
rpSP3.LegendVisible = "on";
rpSP3.LegendLocation = "best";
grid on;
rpSP3.Responses
rpSP3.Responses(1).LineWidth = 2;
rpSP3.Responses(2).LineWidth = 2;
ylim([0, 260]);
rpSP3.AxesStyle(1).FontSize=16;
rpSP3.XLabel(1).FontSize = 16;
rpSP3.YLabel(1).FontSize = 16;
G3info = stepinfo(G3s, "SettlingTimeThreshold", 0.01)
G3approx_info = stepinfo(G3sAPPROX, "SettlingTimeThreshold", 0.01)
Final Example - A Set of Dominant Poles and Zeros
Consider the LTI system described by the transfer function
clear
s = tf('s'); % the transfer function building helper
G4s = 250*(1+0.8*s+(s^2/4))*(1+s/200)/( (1+(s/20))*(1+s/3)*(1+(9/20)*s+s^2/16) ); % the transfer function
set(G4s, "Name", "G4");
G4_staticGAIN = dcgain(G4s) % the static gain of the LTI system
p_G4s = sort(pole(G4s), 'ascend') % the poles of the transfer function
z_G4s = sort(zero(G4s), 'ascend') % the zeros of the transfer function
The poles and zeros configuration of the system is as follows
figure('Units','pixels', 'Position',[1,1,1280, 960]);
hPZP4 = pzplot(G4s);
grid on
% get(hPZP)
% PZPlot properties:
% https://it.mathworks.com/help/ident/ref/controllib.chart.pzplot-properties.html
% https://it.mathworks.com/matlabcentral/answers/2173352-change-marker-colour-of-pzplot
%
hPZP4.Responses(1).MarkerSize = 12;
hPZP4.Responses(1).LineWidth = 2;
hPZP4.Responses(1).Color = 'r';
% get(hPZP.AxesStyle(1))
% see 'doc pzplot' & 'doc pzoptions'
hPZP4.AxesStyle(1).FontSize = 14;
hPZP4.AxesStyle(1).GridColor = [0,0,0];
% hPZP.XLabel(1)
hPZP4.XLabel(1).FontSize = 14;
hPZP4.YLabel(1).FontSize = 14;
The poles closest to the origin are the dominant ones. This time, there are two complex poles and one real pole close to the origin (with real parts very close to each other - we cannot neglect any of the poles in question). In addition, the pair of complex zeros must also be considered.
The Dominant Poles Approximation
We may approximate the system 
using a very simple model: the static gain and the dominant poles
using a very simple model: the static gain and the dominant poles
G4sAPPROX = 250*(1+0.8*s+(s^2/4))/( (1+s/3)*(1+(9/20)*s+s^2/16) ); % the approximating transfer function
set(G4sAPPROX, "Name", "G4_approx");
Let's compare the step response of system with the approximate model:
with the approximate model:hf = figure('Units','pixels', 'Position',[1,1,1280, 960]);
% let's create a blank figure
rpSP4 = stepplot(G4s, G4sAPPROX, [0, 180]);
grid on; zoom on;
rpSP4
% refer to https://it.mathworks.com/help/control/ug/customizing-response-plots-from-the-command-line.html
rpSP4.Characteristics(1).SettlingTime(1).Threshold = 0.01; % 1 %
rpSP4.Characteristics(1).SettlingTime(1).Visible = "on";
rpSP4.Characteristics(1).SteadyState(1).Visible = "on";
rpSP4.LegendVisible = "on";
rpSP4.LegendLocation = "best";
rpSP4.Responses
rpSP4.Responses(1).LineWidth = 2;
rpSP4.Responses(2).LineWidth = 2;
%ylim([0, 260]);
rpSP4.AxesStyle(1).FontSize=16;
rpSP4.XLabel(1).FontSize = 16;
rpSP4.YLabel(1).FontSize = 16;
G4info = stepinfo(G4s, "SettlingTimeThreshold", 0.01)
G4approx_info = stepinfo(G4sAPPROX, "SettlingTimeThreshold", 0.01)