How to fix a code?

1 次查看（过去 30 天）
Valeriy
Valeriy 2025-3-16
移动：Walter Roberson 2025-3-16
I have the implementation code for the optimal speed of an axisymmetric spacecraft turn:
% Оптимальный разворот осесимметричного КА
% 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. Функция для решения системы уравнений
function dxdt = equations(t, x, I1, I3, M3, t_vector)
% x = [q0, q1, q2, q3, omega1, omega2, omega3]
q = x(1:4);
omega = x(5:7);
% Интерполяция управляющего момента
M3_interp = interp1(t_vector, M3, 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
% 6. Численное решение
x0 = [q0; omega0];
options = odeset('RelTol',1e-6,'AbsTol',1e-9);
[T,X] = ode45(@(t,x) equations(t, x, I1, I3, M3, t), 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).
If I try to run it, it will give an error:
Error: File: pere.m Line: 52 Column: 1
Function definitions in a script must appear at the end of the file.
Move all statements after the "equations" function definition to before the first local function definition.
If I move "item 5" to the end of the code, I get this:
Error using interp1>reshapeAndSortXandV
X and V must be of the same length.
Error in interp1 (line 128)
[X,V,orig_size_v] = reshapeAndSortXandV(X,V);
Error in pere>equations (line 81)
M3_interp = interp1(t_vector, M3, t);
Error in pere>@(t,x)equations(t,x,I1,I3,M3,t) (line 31)
[T,X] = ode45(@(t,x) equations(t, x, I1, I3, M3, t), tspan, x0, options);
Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
Error in pere (line 31)
[T,X] = ode45(@(t,x) equations(t, x, I1, I3, M3, t), tspan, x0, options);
Help me fix the code so that it just runs and outputs images.

请先登录，再回答此问题。

回答（2 个）

Star Strider
Star Strider 2025-3-16
Ran in:
I believe it was getting confused by two different uses of ‘t’ in the ode45 argument list. When I added ‘t_vector’ and used that instead, it ran without error. I leave it to you to determine that it is giving the correct result.
Try this —
% Оптимальный разворот осесимметричного КА
% 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. Функция для решения системы уравнений
% 6. Численное решение
x0 = [q0; omega0];
options = odeset('RelTol',1e-6,'AbsTol',1e-9);
t_vector = t;
% figure
% plot(t, M3)
% grid
[T,X] = ode45(@(t,x) equations(t, x, I1, I3, M3, t_vector), 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_final_simulated = 4×1
0.7338 0 0 0.6794
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
q_difference = qf - q_final_simulated;
% 9. Необходима оптимизация времени переключения t_switch для достижения
% минимального времени разворота и точного попадания в конечную
% ориентацию qf и угловую скорость omegaf. Это можно сделать
% используя функции оптимизации MATLAB (например, fminsearch).
function dxdt = equations(t, x, I1, I3, M3, t_vector)
% x = [q0, q1, q2, q3, omega1, omega2, omega3]
q = x(1:4);
omega = x(5:7);
% Интерполяция управляющего момента
M3_interp = interp1(t_vector, M3, 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
.
Walter Roberson
Walter Roberson 2025-3-16
移动：Walter Roberson 2025-3-16
Ran in:
You created a vector named t .
Then in constructing the anonymous function for ode45 purposes, you used
@(t,x) equations(t, x, I1, I3, M3, t)
Notice you used t as your time parameter. But you also used t as your vector, so when you construct the anonymous function, the time parameter t is also used as the time vector t
% Оптимальный разворот осесимметричного КА
% 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. Функция для решения системы уравнений
% 6. Численное решение
x0 = [q0; omega0];
options = odeset('RelTol',1e-6,'AbsTol',1e-9);
[T,X] = ode45(@(T,x) equations(T, x, I1, I3, M3, t), 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_final_simulated = 4×1
0.7338 0 0 0.6794
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
q_difference = qf - q_final_simulated;
% 9. Необходима оптимизация времени переключения t_switch для достижения
% минимального времени разворота и точного попадания в конечную
% ориентацию qf и угловую скорость omegaf. Это можно сделать
% используя функции оптимизации MATLAB (например, fminsearch).
function dxdt = equations(t, x, I1, I3, M3, t_vector)
% x = [q0, q1, q2, q3, omega1, omega2, omega3]
q = x(1:4);
omega = x(5:7);
% Интерполяция управляющего момента
M3_interp = interp1(t_vector, M3, 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

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by