A reference trajectory strategy using PFC

PFC_controller.m

% MATLAB Script for Simulating PFC Controller from Example 4.4

% Clear workspace and command window
clear;
close all;
clc;

% Sampling time
Ts = 3;  % seconds

% Reference trajectory parameter
alpha = exp(-1/3);  % Approximately 0.7165

% Simulation parameters
N = 15;  % Total simulation steps
r = 3;  % Set-point (constant at 3)

% Initialize variables
y = zeros(N+1, 1);     % Plant output
u = zeros(N+1, 1);     % Control input
w = zeros(N+1, 2);     % Reference trajectory
yf = zeros(N+1, 2);    % Free response
Delta_u = zeros(N+1, 1); % Control increments

% Initial conditions
yf(1,1) = 2;
yf(1,2) = 2;
u(1) = 0.3;
% Diophantine equation coefficients (from Example 4.4)
f10 = 0.7;
f20 = 0.49;

% Step responses
g1 = 2;           % g(1)
g2 = 3.4;         % g(2)
G = [g1; g2];               % G matrix

y(1) = 2;

% Main simulation loop
for t = 1:N+1
    % Compute the reference trajectory at t+1 and t+2
    w(t, 1) = r - alpha * (r - y(t));
    w(t, 2) = r - alpha^2 * (r - y(t));

    T = [w(t, 1); w(t, 2)];

    % Compute free responses
    % For k = 1
    yf(t,1) = f10 * y(t) + g1 * u(t);
    % For k = 2
    yf(t,2) = f20 * y(t) + g2 * u(t);

    % Construct (T - Yf) vector
    T_minus_Yf = [w(t, 1) - yf(t,1); w(t,2) - yf(t,2)];

    % Solve for control increment (Delta u)
    Delta_u(t) = (G' * G) \ (G' * T_minus_Yf);

    % Update control input
    u(t+1) = u(t) + Delta_u(t);
    % Update plant output
    y(t+1) = 0.7 * y(t) + 2 * u(t+1);
end


% Time vector for plotting
time = (0:N+1) * Ts;

% Extract data for plotting
y_plot = y(1:N+2);
u_plot = u(1:N+2);
time_plot = time;

% Prepare w for plotting
% Since w is computed for t = 1:N+1, we'll create a time vector for w
w_time = time(1:end-1) + Ts;  % w(t) corresponds to time t + 1

% Plotting the results
figure;

% Plot Plant Output vs. Set-Point and Reference Trajectory
subplot(2,1,1);
plot(time_plot, y_plot, 'b-o', 'LineWidth', 2);
hold on;
plot(time_plot, r * ones(size(time_plot)), 'r--', 'LineWidth', 1.5);
plot(w_time, w(:,1), 'g-*', 'LineWidth', 1.5);  % Plot w(t+1)
xlabel('Time (seconds)');
ylabel('Output y(t)');
title('Plant Output, Set-Point, and Reference Trajectory');
legend('Plant Output y(t)', 'Set-Point r(t)', 'Reference Trajectory w(t+1)', 'Location', 'Best');
grid on;

% Plot Control Input over Time
subplot(2,1,2);
stairs(time_plot, u_plot, 'k-o', 'LineWidth', 2);
xlabel('Time (seconds)');
ylabel('Control Input u(t)');
title('Control Input over Time');
grid on;