hi.I have problem in spacecraft dynamic control simulation.

someone PLS help me the part of my code is this

global Ix Kx Kxd phi_com Iy Ky Kyd theta_com Kz Kzd psi_com Iz w0;
Ix          = 9;
Kx          = 9;
Kxd         = -42;
phi_com     = 15 * pi / 180;
Iy          = 13;
Ky          = 13;
Kyd         = -91;
theta_com   = 15 * pi / 180;
Iz          = 13;
Kz          = 13;
Kzd         = -91; 
psi_com     = 15 * pi / 180;
tmax        = 50;
y0          = [0; 0; 0; 0; 0; 0; 0];
w0          = 0.001038;

it runs well when the kxd and kyd and kzd are positive but when I negative them as U can see, the code doesnt run I dont know what happens I have to negative them because of my refrence article. other parts of the code are correct and run well when they have positive amount untill I negative them

[Sam Chak]

@ali.heydari, The code displayed only involves assigning values to the variables. However, the sections where you define the functions for (1) Spacecraft Dynamics and (2) the Attitude Control Law are not included. It would be helpful if you could provide that two separate functions as well.

function dydt = SpacecraftDynamics(t, y)
    ...
end

function tau  = AttitudeController(t, y)
    ...
end

[ali.heydari]

let me send U all the code which is made of 4 parts

global Ix Kx Kxd phi_com Iy Ky Kyd theta_com Kz Kzd psi_com Iz w0;

Ix = 9;

Kx = 9;

Kxd = -42;

phi_com = 15 * pi / 180;

Iy = 13;

Ky = 13;

Kyd = -91;

theta_com = 15 * pi / 180;

Iz = 13;

Kz = 13;

Kzd = -91;

psi_com = 15 * pi / 180;

tmax = -50;

y0 = [0; 0; 0; 0; 0; 0; 0];

w0 = 0.001038;

first part

function xdot = ode(~, y)

global Ix Kx Kxd phi_com Iy Ky Kyd theta_com Kz Kzd psi_com Iz w0

phi = y(1);

theta = y(2);

psi = y(3);

wx = y(4);

wy = y(5);

wz = y(6);

AS = DCM(phi, theta, psi); % Spacecraft Direction Cosine Matrix

AT = DCM(phi_com, theta_com, psi_com); % Target DCM

qt4 = 0.5 * sqrt(1 + AT(1,1) + AT(2,2) + AT(3,3)); % Page 325 Eq A.4.16 (Appendix)

qt1 = (0.25 * (AT(2,3) - AT(3,2))) / qt4;

qt2 = (0.25 * (AT(3,1) - AT(1,3))) / qt4;

qt3 = (0.25 * (AT(1,2) - AT(2,1))) / qt4;

qT = [qt4 qt3 -qt2 qt1;

-qt3 qt4 qt1 qt2;

qt2 -qt1 qt4 qt3;

-qt1 -qt2 -qt3 qt4];

qs4 = 0.5 * sqrt(1 + AS(1,1) + AS(2,2) + AS(3,3)); % Page 325 Eq A.4.16 (Appendix)

qs1 = (0.25 * (AS(2,3) - AS(3,2))) / qs4;

qs2 = (0.25 * (AS(3,1) - AS(1,3))) / qs4;

qs3 = (0.25 * (AS(1,2) - AS(2,1))) / qs4;

qS = [-qs1; -qs2; -qs3; qs4];

qE = qT * qS; % Quaternion Error Vector

q1E = qE(1); % Extracting quaternion error components

q2E = qE(2);

q3E = qE(3);

q4E = qE(4);

WRI = [0; -w0; 0];

WRIB = AS * WRI;

W = [wx; wy; wz];

WBR = W - WRIB;

p = WBR(1);

q = WBR(2);

r = WBR(3);

phidot = p + (q * sin(phi) + r * cos(phi)) * tan(theta); % 3-2-1

thetadot = q * cos(phi) - r * sin(phi);

psidot = ((q * sin(phi)) + r * cos(phi)) * sec(theta);

Tcx = (2 * Kx) * (q1E * q4E) - (Kxd) * p; % Page 156 Eqs 7.2.15

Tcy = (2 * Ky) * (q2E * q4E) - (Kyd) * q;

Tcz = (2 * Kz) * (q3E * q4E) - (Kzd) * r;

wxdot = (Tcx - wy * wz * (Iz - Iy)) / (Ix); % Page 95 Eq 4.5.1 Euler's Moment Equations with a little bit Modification

wydot = (Tcy - wx * wz * (Ix - Iz)) / (Iy);

wzdot = (Tcz - wy * wx * (Iy - Ix)) / (Iz);

AE = AS * AT';

alpha = acos(0.5 * (trace(AE) - 1));

Eulerintdot = alpha * 180 / pi; % Page 157

xdot = [phidot; thetadot; psidot; wxdot; wydot; wzdot; Eulerintdot];

end

second part

function AS = DCM(phi, theta, psi)

Aphi = [1 0 0;

0 cos(phi) sin(phi);

0 -sin(phi) cos(phi)]; % Page 321 Eq A.3.4

Atheta = [cos(theta) 0 -sin(theta);

0 1 0;

sin(theta) 0 cos(theta)]; % Page 321 Eq A.3.3

Apsi = [cos(psi) sin(psi) 0;

-sin(psi) cos(psi) 0;

0 0 1]; % Page 321 Eq A.3.2

AS = Aphi * Atheta * Apsi;

end

third part

clc

clear;

close all;

initial

tmax=50;

[t, y] = ode45(@ode, [0 tmax], y0);

% Plot results

figure;

subplot(3,1,1);

plot(t, y(:,4), 'linewidth', 2); % Plot P

xlabel('Time (sec)');

ylabel('P');

title('Angular Velocities');

subplot(3,1,2);

plot(t, y(:,5), 'linewidth', 2); % Plot Q

xlabel('Time (sec)');

ylabel('Q');

subplot(3,1,3);

plot(t, y(:,6), 'linewidth', 2); % Plot R

xlabel('Time (sec)');

ylabel('R');

figure;

subplot(3,1,1);

plot(t, y(:,1)*180/pi, 'linewidth', 2); % Plot phi

xlabel('Time (sec)');

ylabel('Phi');

title('Euler Angles');

subplot(3,1,2);

plot(t, y(:,2)*180/pi, 'linewidth', 2); % Plot theta

xlabel('Time (sec)');

ylabel('Theta');

subplot(3,1,3);

plot(t, y(:,3)*180/pi, 'linewidth', 2); % Plot psi

xlabel('Time (sec)');

ylabel('Psi');

fourth part

it doesnt run

I dont know if understand your words well or not

Answer 1

[Sam Chak]

Hi @ali.heydari

I renamed the ODE function to 'odefcn()' because 'ode' is a built-in function introduced in the R2023b release. In your original code, there were two tmax parameters. By disabling the one with a negative value, the code is now running properly. I have simulated the system up to 2.5 seconds to provide you with the results, but I didn't 'beautify' the code.

Update: I have made adjustments to the control algorithm based on the equations provided in Sidi's book. As a result, the attitude control system is now stable.

global Ix Kx Kxd phi_com Iy Ky Kyd theta_com Kz Kzd psi_com Iz w0;
Ix = 9;
Kx = 9;
Kxd = -42;
phi_com  = 15 * pi / 180;
Iy = 13;
Ky = 13;
Kyd = -91;
theta_com = 15 * pi / 180;
Iz = 13;
Kz = 13;
Kzd = -91; 
psi_com   = 15 * pi / 180;
% tmax = -50;
y0 = [0; 0; 0; 0; 0; 0; 0];
w0 = 0.001038;
    
function xdot = odefcn(~, y)
    global Ix Kx Kxd phi_com Iy Ky Kyd theta_com Kz Kzd psi_com Iz w0
    phi = y(1);
    theta = y(2);
    psi = y(3);
    wx = y(4);
    wy = y(5);
    wz = y(6);
    AS = DCM(phi, theta, psi); % Spacecraft Direction Cosine Matrix
    AT = DCM(phi_com, theta_com, psi_com); % Target DCM
    qt4 = 0.5 * sqrt(1 + AT(1,1) + AT(2,2) + AT(3,3)); % Page 325 Eq A.4.16 (Appendix)
    qt1 = (0.25 * (AT(2,3) - AT(3,2))) / qt4;
    qt2 = (0.25 * (AT(3,1) - AT(1,3))) / qt4;
    qt3 = (0.25 * (AT(1,2) - AT(2,1))) / qt4;
    qT = [qt4 qt3 -qt2 qt1;
          -qt3 qt4 qt1 qt2;
          qt2 -qt1 qt4 qt3;
          -qt1 -qt2 -qt3 qt4];
    qs4 = 0.5 * sqrt(1 + AS(1,1) + AS(2,2) + AS(3,3)); % Page 325 Eq A.4.16 (Appendix)
    qs1 = (0.25 * (AS(2,3) - AS(3,2))) / qs4;
    qs2 = (0.25 * (AS(3,1) - AS(1,3))) / qs4;
    qs3 = (0.25 * (AS(1,2) - AS(2,1))) / qs4;
    qS = [-qs1; -qs2; -qs3; qs4];
    qE = qT * qS; % Quaternion Error Vector
    q1E = qE(1); % Extracting quaternion error components
    q2E = qE(2);
    q3E = qE(3);
    q4E = qE(4);
    WRI = [0; -w0; 0];
    WRIB = AS * WRI;
    W = [wx; wy; wz];
    WBR = W - WRIB;
    p = WBR(1);
    q = WBR(2);
    r = WBR(3);
    phidot = p + (q * sin(phi) + r * cos(phi)) * tan(theta); % 3-2-1
    thetadot = q * cos(phi) - r * sin(phi);
    psidot = ((q * sin(phi)) + r * cos(phi)) * sec(theta);
    
    %% ------ changes made ------
    Tcx = (2 * Kx) * (q1E * q4E) + (Kxd) * p; % Page 156 Eqs 7.2.15
    Tcy = (2 * Ky) * (q2E * q4E) + (Kyd) * q;
    Tcz = (2 * Kz) * (q3E * q4E) + (Kzd) * r;
    %% ------ changes made ------
    
    wxdot = (Tcx - wy * wz * (Iz - Iy)) / (Ix); % Page 95 Eq 4.5.1 Euler's Moment Equations with a little bit Modification
    wydot = (Tcy - wx * wz * (Ix - Iz)) / (Iy);
    wzdot = (Tcz - wy * wx * (Iy - Ix)) / (Iz);
    AE = AS * AT';
    alpha = acos(0.5 * (trace(AE) - 1));
    Eulerintdot = alpha * 180 / pi; % Page 157
    xdot = [phidot; thetadot; psidot; wxdot; wydot; wzdot; Eulerintdot];
end

function AS = DCM(phi, theta, psi)
    Aphi = [1 0 0;
            0 cos(phi) sin(phi);
            0 -sin(phi) cos(phi)]; % Page 321 Eq A.3.4
    Atheta = [cos(theta) 0 -sin(theta);
              0 1 0;
              sin(theta) 0 cos(theta)]; % Page 321 Eq A.3.3
    Apsi = [cos(psi) sin(psi) 0;
            -sin(psi) cos(psi) 0;
            0 0 1]; % Page 321 Eq A.3.2
    AS = Aphi * Atheta * Apsi;
end   

tmax = 50;
[t, y] = ode45(@odefcn, [0 tmax], y0);

% Plot results
figure;
subplot(3,1,1);
plot(t, y(:,4), 'linewidth', 2); % Plot P
xlabel('Time (sec)');
ylabel('P');
title('Angular Velocities');
subplot(3,1,2);
plot(t, y(:,5), 'linewidth', 2); % Plot Q
xlabel('Time (sec)');
ylabel('Q');
subplot(3,1,3);
plot(t, y(:,6), 'linewidth', 2); % Plot R
xlabel('Time (sec)');
ylabel('R');

figure;
subplot(3,1,1);
plot(t, y(:,1)*180/pi, 'linewidth', 2); % Plot phi
xlabel('Time (sec)');
ylabel('Phi');
title('Euler Angles');
subplot(3,1,2);
plot(t, y(:,2)*180/pi, 'linewidth', 2); % Plot theta
xlabel('Time (sec)');
ylabel('Theta');
subplot(3,1,3);
plot(t, y(:,3)*180/pi, 'linewidth', 2); % Plot psi
xlabel('Time (sec)');
ylabel('Psi');

[ali.heydari]

thanks my friend love U so much but in the refrence that I use the plots are reverse let me send u the picture of the plot. what can I do to reverse the plots??

[ali.heydari]

sorry if I answere the next message late here it is 3 A.M I worked hard and feel asleep.see U tomorrow