function Principal

% Example code showing how to compute and plot the principal gains of a system.

% Define the system matrices.
A = [-1 0; 0 -1];
B = [1 0; 0 1];
C = [1 2; 2 1];
D = [0 0; 0 0];
I = eye(size(A));
% Define the frequency variable.
w = logspace(-2, 2, 100);
% Compute the principal gains.
for k = 1 : 100
    G = C * (j*w(k)*I - A)^(-1) * B + D;
    Pgains(:, k) = svd(G);
end
% Plot the principal gains.
close all;
semilogx(w, 20*log10(Pgains));
xlabel('Frequency (rad/s)');
ylabel('Principal gains (dB)');

% An easier way to plot the principal gains.
figure;
sigma(A, B, C, D);

% Get the principal gains at 10 rad/s.
sys = ss(A, B, C, D);
G = evalfr(sys, 10i);
[U, S, V] = svd(G);
disp(['Principal gains at 10 rad/s = ', num2str(S(1,1)), ', ', num2str(S(2,2))]);  

% Simulation
t = 0 : 0.01 : 10;
u(1,:) = 1*sin(10*t); % arbitrary input
u(2,:) = 4*sin(10*t-10);

%u(1,:) = V(1,1)*sin(10*t); % maximizing input
%u(2,:) = V(2,1)*sin(10*t);

%u(1,:) = V(1,2)*sin(10*t); % minimizing input
%u(2,:) = V(2,2)*sin(10*t);

[y, t] = lsim(sys, u, t);
figure;
subplot(2,1,1);
plot(t, u(1,:), t, u(2,:));
ylabel('input');
subplot(2,1,2);
plot(t, y(:,1), t, y(:,2));
ylabel('output');
xlabel('time');