clear all;
clc;
% Parámetros del modelo
m = 5.8; % 58000 toneladas = 5.8e7 kg
b = 2.5; % 25000 Ns/m = 2.5e4 Ns/m
% Parámetros de la ley adaptativa
a2 = m;
a1 = b;
a0 = 0; % Ganancia
% Parámetros del modelo de referencia
Ts = 0.0556; % Periodo natural (horas)
zet = 1; % Tasa de amortiguamiento
wn = 2 * pi / Ts; % Frecuencia natural
am2 = 1 / wn^2;
am1 = 2 * wn * zet * am2;
am0 = 1;
% Parámetros del controlador
we = 0.5 * wn; % Mitad de la frecuencia natural
zet_e = 0.7; % Amortiguamiento del controlador
bet1 = 2 * we * zet_e;
bet0 = we^2;
% Parámetros de la ley adaptativa
gam = [1e-3 0; 0 1];
Q = 0.113 * eye(2);
A = [0 1; -bet0 -bet1];
P = lyap(A, Q);
P11 = P(1,1);
P22 = P(2,2);
% Simulación
t = 0:0.1:12; % Tiempo en horas
r = 20 * square(2 * pi * 1 / 2 * t); % Referencia cuadrada
y = zeros(size(t)); % Salida inicial
ym = zeros(size(t)); % Modelo de referencia
u = zeros(size(t)); % Entrada de control
e = zeros(size(t)); % Error
% Inicialización de valores iniciales
e(1) = r(1); % Error inicial
ym(1) = 0; % Salida inicial del modelo de referencia
y(1) = 0; % Salida inicial del sistema
for k = 2:length(t)
e(k) = r(k) - y(k-1); % Error
ym(k) = am2 * e(k-1); % Modelo de referencia (simplificado para evitar índices negativos)
u(k) = -P11 * e(k-1) - P22 * e(k-1); % Control adaptativo
y(k) = a2 * u(k) + a1 * e(k); % Salida del sistema
end
% Gráficos
figure;
subplot(2, 2, 1);
plot(t, r, '--', t, y, 'LineWidth', 1.5);
xlabel('Time (hours)');
ylabel('Tracking (meters)');
legend('Reference', 'Output');
title('Tracking Performance');
subplot(2, 2, 2);
plot(t, e, 'LineWidth', 1.5);
xlabel('Time (hours)');
ylabel('Error');
title('Error Signal');
subplot(2, 2, 3);
plot(t, u, 'LineWidth', 1.5);
xlabel('Time (hours)');
ylabel('Control Effort');
title('Control Effort');
subplot(2, 2, 4);
plot(t, ym, 'LineWidth', 1.5);
xlabel('Time (hours)');
ylabel('Model Output');
title('Model Reference');