Исправление программы MATLAB

 % Оптимальный разворот осесимметричного КА  
 
% 1. Параметры КА  
I1 = 100;   % Момент инерции относительно оси X и Y (kg*m^2)  
I3 = 50;    % Момент инерции относительно оси Z (kg*m^2)  
M3_max = 1; % Максимальный управляющий момент (Nm)  
 
% 2. Начальные и конечные условия  
q0 = [1; 0; 0; 0]; % Начальный кватернион (ориентация)  
omega0 = [0; 0; 0.1]; % Начальная угловая скорость (rad/s)  
qf = [0.7071; 0; 0; 0.7071]; % Конечный кватернион (разворот на 90 градусов вокруг X)  
omegaf = [0; 0; 0];  % Конечная угловая скорость  
 
% 3. Параметры моделирования  
tspan = [0 20]; % Временной интервал (секунды)  
dt = 0.01;      % Шаг по времени  
t = tspan(1):dt:tspan(2);  
 
% 4.  Оптимальное управление (Пример реализации с переключением по времени,  
%     необходимо оптимизировать моменты переключения!)  
t_switch = 5; % Пример времени переключения  
M3 = zeros(size(t));  
M3(t < t_switch) = M3_max;  
M3(t >= t_switch) = -M3_max;  
 
% 5. Анонимная функция для системы уравнений 
odefun = @(t, x) equations_func(t, x, I1, I3, M3, t); 
 
% 6. Численное решение  
x0 = [q0; omega0];  
options = odeset('RelTol',1e-6,'AbsTol',1e-9);  
[T,X] = ode45(odefun, tspan, x0, options);  
 
% 7. Извлечение результатов  
Q = X(:,1:4); % Кватернионы  
Omega = X(:,5:7); % Угловые скорости  
 
% 8.  Визуализация результатов  
 
% 8.1  Угловая скорость Omega_3  
figure;  
plot(T, Omega(:,3));  
xlabel('Время (с)');  
ylabel('Угловая скорость Omega_3 (рад/с)');  
title('Изменение угловой скорости Omega_3');  
grid on;  
 
% 8.2  Фазовый портрет Omega_3 - Omega_3' (приблизительно)  
omega3_dot = diff(Omega(:,3))/dt; % Численное дифференцирование  
figure;  
plot(Omega(1:end-1,3), omega3_dot);  
xlabel('Omega_3 (рад/с)');  
ylabel('Omega_3'' (рад/с^2)');  
title('Фазовый портрет Omega_3');  
grid on;  
 
% 8.3  Трёхмерная траектория угловой скорости (фазовое пространство)  
figure;  
plot3(Omega(:,1), Omega(:,2), Omega(:,3));  
xlabel('Omega_1 (рад/с)');  
ylabel('Omega_2 (рад/с)');  
zlabel('Omega_3 (рад/с)');  
title('Траектория угловой скорости в фазовом пространстве');  
grid on;  
 
% 8.4  Проверка достижения конечной ориентации (примерно)  
q_final_simulated = Q(end,:)';  
q_difference = qf - q_final_simulated;  
 
% 9.  Необходима оптимизация времени переключения t_switch для достижения  
%     минимального времени разворота и точного попадания в конечную  
%     ориентацию qf и угловую скорость omegaf.  Это можно сделать  
%     используя функции оптимизации MATLAB (например, fminsearch). 
 
% Вспомогательная функция для уравнений движения 
function dxdt = equations_func(t, x, I1, I3, M3_values, t_values) 
    % x = [q0, q1, q2, q3, omega1, omega2, omega3]  
    q = x(1:4);  
    omega = x(5:7);  
 
    % Интерполяция управляющего момента  
    M3_interp = interp1(t_values, M3_values, t); 
 
    % Уравнения Эйлера  
    omega_dot = zeros(3,1);  
    omega_dot(1) = ((I1 - I3)/I1) * omega(2) * omega(3);  
    omega_dot(2) = ((I3 - I1)/I1) * omega(3) * omega(1);  
    omega_dot(3) = M3_interp / I3;  
 
    % Уравнение для кватерниона  
    Omega_matrix = [  0    -omega(1)   -omega(2)   -omega(3);  
                     omega(1)    0     omega(3)   -omega(2);  
                     omega(2)   -omega(3)    0     omega(1);  
                     omega(3)    omega(2)   -omega(1)    0    ];  
    q_dot = 0.5 * Omega_matrix * q;  
 
    dxdt = [q_dot; omega_dot];  
end