Robustness
Housekeeping
close all; clear; clc; s = tf('s');
Suppose we wish to analyze a mass-spring-damper system with:
x = position
m = mass
b = damper constant
k = spring constant
u = input force
The transfer function then will be: 
However, consider a realistic situation where the values are not single values, but rather given as normal distributions with a mean and standard deviation.
t = 0:0.1:10;
plant_zeros = [];
plant_poles = [];
yout = zeros(size(t));
i=1;
for i = 1:1:100
m = 0.2*randn(1,1)+10;
b = 0.2*randn(1,1)+5;
k = 0.2*randn(1,1)+3;
Plant = (1)/(m*s*s + b*s + k);
Controller = 1;
System = feedback(Controller*Plant,1);
plant_zeros = [plant_zeros zero(System)];
plant_poles = [plant_poles pole(System)];
[y, t] = step(System, t);
yout = [yout; y'];
end
f = figure();
f.Position = [0 0 1500 500];
sgtitle('Step Response Based on Pole');
subplot(1,2,1)
title('Pole-Zero Map');
scatter(real(plant_poles), imag(plant_poles), 'Marker', 'x', 'MarkerEdgeColor','k');
xlabel('frequency (rad/s)');
ylabel('magnitude');
subplot(1,2,2)
for i = 1:1:100
plot(t, yout(i,:)); hold on;
end
xlabel('frequency (rad/s)');
ylabel('phase (rad)');
hold off;
t = 0:0.1:10;
plant2_zeros = [];
plant2_poles = [];
yout2 = zeros(size(t));
i=1;
for i = 1:1:100
m2 = 0.2*randn(1,1)+10;
b2 = 0.2*randn(1,1)+5;
k2 = 0.2*randn(1,1)+3;
Plant2 = (1)/(m2*s*s + b2*s + k2);
Controller2 = 50; % SET THE CONTROLLER HERE!
System2 = feedback(Controller2*Plant2,1);
plant2_zeros = [plant2_zeros zero(System2)];
plant2_poles = [plant2_poles pole(System2)];
[y2, t] = step(System2, t);
yout2 = [yout2; y2'];
end
f = figure();
f.Position = [0 0 1500 500];
sgtitle('Step Response Based on Pole');
subplot(1,2,1)
title('Pole-Zero Map');
scatter(real(plant_poles), imag(plant_poles), 'Marker', 'x', 'MarkerEdgeColor','r'); hold on;
scatter(real(plant2_poles), imag(plant2_poles), 'Marker', 'x', 'MarkerEdgeColor','g'); hold off;
xlabel('frequency (rad/s)');
ylabel('magnitude');
legend('red=original', 'green=new', 'FontSize', 10, 'Location', 'northeast')
subplot(1,2,2)
for i = 1:1:100
plot(t, yout2(i,:)); hold on;
end
xlabel('frequency (rad/s)');
ylabel('phase (rad)');
hold off;
