function PendRGA

% Compute the relative gain array elements of the inverted pendulum system.
% The system inputs are horizontal force and pendulum acceleration.
% The system outputs are horizontal distance and pendulum angle.

m = 0.2; % pendulum mass
M = 1.0; % cart mass
g = 9.81; % acceleration due to gravity
F = 0.1; % coefficient of viscous friction on cart wheels
l = 1; % length of pendulum
r = 0.02; % radius of pendulum mass
J = m * r * r / 2; % moment of inertia of pendulum mass (cylinder)

% Linearized system matrices.
A = [0       1             0   0;
     0    -F/M        -m*g/M   0;
     0       0             0   1;
     0   F/M/l   (M+m)*g/M/l   0];
B = [0 0 ; 1/M 0 ; 0 0 ; -1/M/l 1];
C = [1 0 0 0 ; 0 0 1 0];
D = [0 0 ; 0 0];

% Compute the RGA.
sys = ss(A, B, C, D);
nw = 200;
w = logspace(-2, 2, nw);
h = freqresp(sys, w);
for i = 1 : nw
    hinv = pinv(h(:,:,i)).';
    RGA(:,:,i) = h(:,:,i) .* hinv;
end
% Plot the RGA elements.
close all;
semilogx(w, squeeze(RGA(1,1,:)), 'r:', w, squeeze(RGA(1,2,:)), 'b', w, squeeze(RGA(2,1,:)), 'b', w, squeeze(RGA(2,2,:)), 'r:');
xlabel('frequency (rad/s)');
ylabel('interactions');
legend('RGA(1,1) = RGA(2,2)', 'RGA(1,2) = RGA(2,1)');