ARGOS

Problem Statement

Current satellites are not reusable, and recovery is costly. As end-of-life satellites accumulate, space debris poses increasing risks to active missions. Developing identifiable, repairable, and serviceable spacecraft requires new research, but tools to replicate the tumble, rotation, and motion of unresponsive satellites are limited. ARGOS provides a foundational test platform to simulate these conditions, enabling research into satellite reusability and debris mitigation strategies.

Bill of Materials (BOM)

The following table lists the components used in the prototype, including part numbers, quantities, materials, estimated costs, and potential suppliers.

Stewart Platform Evaluation MATLAB Script

PYTHON
close all
clear all
clc

%======================= INITIAL PARAMETERS =======================
%----- Actuator Properties -----
linear_vel = 0.04;                   % Linear speed of actuator (m/s)
axial_max = 150;                     % Max axial force allowed (N)
U_joint = 0.08;                      % Height of U-joint (m)
L0 = 0.290 + 2*U_joint;              % Retracted length (m)
stroke = 0.150;                      % Stroke length (m)
act_d = (0.787)*0.0254;              % Actuator stroke bar diameter (in → m)
L_max = L0 + stroke;                 % Max actuator length (m)
E = 68e9;                            % Modulus of elasticity (Pa)
YS = 200e6;                          % Yield strength (Pa)
rot_act_h = [0 0 0.0395]';           % Rotary actuator height (m)

%----- Plate Geometry -----
r_top_act = 0.1;                     % Top plate actuator radius (m)
r_bot_act = 0.21;                    % Bottom plate actuator radius (m)
space = 0.03;                        % Distance from actuator to plate edge (m)
twist_clock = -15;                   % Rotation between top and bottom plates (deg)
twist_counter = 15;                  % Rotation between top and bottom plates (deg)
rho_top = 2700;                      % Density of top plate material (kg/m³)
V_top = 4.25e-3;                     % Volume of top plate (m³)
rho_bot = 2700;                      % Density of bottom plate material (kg/m³)
V_bot = 5e-3;                        % Volume of bottom plate (m³)

%----- Satellite Bus -----
a = 0.460; b = 0.460; c = 0.460;     % Dimensions of bus (m)
rho_bus = 1240;                      % Density of bus material (kg/m³)
t = 0.050;                           % Wall thickness (m)
cushion = 0.005;                     % Clearance for contact (m)

%======================= ROTATION SETUP ============================
roll_range  = 0:1:45;                % Roll angles (deg)
pitch_range = 0:1:45;                % Pitch angles (deg)
yaw_range   = 0:1:360;               % Yaw angles (deg)

%======================= STORAGE INITIALIZATION ===================
SafeDesigns = [];
SafeFOS = [];
SafeOmega = [];
FailCases_all = {};
safeCombos = [];
safeCycles = [];

%======================= GEOMETRY SWEEP CONTROL ===================
violation_flag = false;
failureFlag = false;
FailCases = [];
cycleComplete = true;
idx = 1;

% Compute angular speed for given geometry
omega = (linear_vel / r_top_act) * (180/pi); % deg/s

% Plate actuator locations (bottom plate) (deg)
% Clockwise bottom plate angles
act_theta1_bot = 60;
act_theta2_bot = 180;
act_theta3_bot = 300;
% Counter-Clockwise bottom plate angles
act_theta4_bot = 0;
act_theta5_bot = 120;
act_theta6_bot = 240;

% Plate position setup (clockwise)
[a1_body_bot,a2_body_bot,a3_body_bot,a1_body_top,a2_body_top,a3_body_top] = ...
    PlatePos(act_theta1_bot,act_theta2_bot,act_theta3_bot,r_top_act,r_bot_act,twist_clock);

% Plate position setup (counter-clockwise) -- FIXED to correct indices
[a4_body_bot,a5_body_bot,a6_body_bot,a4_body_top,a5_body_top,a6_body_top] = ...
    PlatePos(act_theta4_bot,act_theta5_bot,act_theta6_bot,r_top_act,r_bot_act,twist_counter);

%---------------------- Weight Calculations ------------------------
W_top = 9.81 * rho_top * V_top;
W_bot = 9.81 * rho_bot * V_bot;
W_sat = 110;                          % Satellite weight (N)
W_rot = (2.5) * 9.81;                 % Rotary actuator weight (N)
W_total = W_top + W_sat + W_rot;      % Total weight on actuators (N)

%---------------------- Result Storage Structure ------------------
Results = struct('maxRoll', [], 'MaxPitch', [], ...
    'roll', [], 'pitch', [], 'yaw', [], ...
    'FOS_yield', [], 'FOS_buckling', [], 'FOS_axial', [], ...
    'dist', [], 'pass', [], 'F1', [], 'F2', [], 'F3', [], ...
    'F4', [], 'F5', [], 'F6', [], 'Rot_act_global', [], ...
    'Bus_CM_global', [], 'L', [], 'Uvec', []);   % <-- added L and Uvec

%======================= ROTATION SWEEP ============================
for phi = roll_range
    for theta_DCM = pitch_range
        % This flag means "did any yaw fail?" — we require ALL yaw to pass
        anyYawFailed = false;

        % Precompute platform rotation matrix (roll & pitch)
        Rx = [1 0 0; 0 cosd(phi) -sind(phi); 0 sind(phi) cosd(phi)];
        Ry = [cosd(theta_DCM) 0 sind(theta_DCM); 0 1 0; -sind(theta_DCM) 0 cosd(theta_DCM)];
        R_platform = Ry * Rx;

        % Compute feasible h0 range using plate geometry only (plate geometry doesn't depend on yaw)
        dxy_all = [ ...
            sqrt((a1_body_top(1)-a1_body_bot(1))^2 + (a1_body_top(2)-a1_body_bot(2))^2), ...
            sqrt((a2_body_top(1)-a2_body_bot(1))^2 + (a2_body_top(2)-a2_body_bot(2))^2), ...
            sqrt((a3_body_top(1)-a3_body_bot(1))^2 + (a3_body_top(2)-a3_body_bot(2))^2), ...
            sqrt((a4_body_top(1)-a4_body_bot(1))^2 + (a4_body_top(2)-a4_body_bot(2))^2), ...
            sqrt((a5_body_top(1)-a5_body_bot(1))^2 + (a5_body_top(2)-a5_body_bot(2))^2), ...
            sqrt((a6_body_top(1)-a6_body_bot(1))^2 + (a6_body_top(2)-a6_body_bot(2))^2)];

            h0_lower_each = sqrt(max(0, L0^2 - dxy_all.^2));
            h0_upper_each = sqrt(max(0, L_max^2 - dxy_all.^2));
            h0_min_feasible = max(h0_lower_each);
            h0_max_feasible = min(h0_upper_each);


        if h0_min_feasible > h0_max_feasible
            % Geometry infeasible for this (phi,theta)
            FailCases(end+1,:) = [-1, phi, theta_DCM, NaN, NaN, NaN];
            continue
        end

        h0 = 0.5 * (h0_min_feasible + h0_max_feasible);
        top_pos = [0 0 h0]';

        % Compute top plate actuator global positions / actuator vectors (independent of yaw)
        [a1_global_top,a2_global_top,a3_global_top,v1,v2,v3,L1,u1,u2,u3,angle1,angle2,angle3] = ...
            PlateRot(R_platform,top_pos,a1_body_top,a2_body_top,a3_body_top,a1_body_bot,a2_body_bot,a3_body_bot);

        [a4_global_top,a5_global_top,a6_global_top,v4,v5,v6,L2,u4,u5,u6,angle4,angle5,angle6] = ...
            PlateRot(R_platform,top_pos,a4_body_top,a5_body_top,a6_body_top,a4_body_bot,a5_body_bot,a6_body_bot);

        % Stroke limit check (all 6 actuators)
        L_all = [L1(:); L2(:)];
        if any((L_all < L0) | (L_all > L_max))
            % stroke infeasible
            continue
        end

        % Precompute constant pieces
        CM_rot_act = 0.5 * rot_act_h;
        Bus_pos0 = [0 0 c/2]';

        % Loop through all yaw values; require ALL yaw to pass
        for psi = yaw_range
            Rz = [cosd(psi) -sind(psi) 0; sind(psi) cosd(psi) 0; 0 0 1];
            R_yaw = Rz * R_platform; % yaw applied on top of platform orientation

            % Rot_act_global and Bus_CM_global rotate with yaw
            Rot_act_global = top_pos + R_yaw * CM_rot_act;
            Bus_CM_global = top_pos + R_yaw * (rot_act_h + Bus_pos0);

            %=== Static equilibrium for all 6 actuators at once ===
            [m1,m2,m3,m4,m5,m6, A, Mom_sat, Mom_rot,F1,F2,F3,F4,F5,F6, x_reac, y_reac, M_tot, M_act] = ...
                Stability6(a1_body_bot,a2_body_bot,a3_body_bot, a4_body_bot,a5_body_bot,a6_body_bot, ...
                           u1,u2,u3,u4,u5,u6, Bus_CM_global,W_sat,W_rot,Rot_act_global,W_total,W_bot);

            %=== Structural analysis (use the triple-wise Analysis as you have) ===
            [sigma_bend_1,sigma_bend_2,sigma_bend_3,deflect_1,deflect_2,deflect_3, ...
             tau_shear_1,tau_shear_2,tau_shear_3,sigma_comb_1,sigma_comb_2,sigma_comb_3, ...
             Pcr_1,Pcr_2,Pcr_3] = Analysis(F1,F2,F3,L1,act_d,E,u1,u2,u3,W_total);

            [sigma_bend_4,sigma_bend_5,sigma_bend_6,deflect_4,deflect_5,deflect_6, ...
             tau_shear_4,tau_shear_5,tau_shear_6,sigma_comb_4,sigma_comb_5,sigma_comb_6, ...
             Pcr_4,Pcr_5,Pcr_6] = Analysis(F4,F5,F6,L2,act_d,E,u4,u5,u6,W_total);

            % Buckling FOS (6 actuators)
            FOS_buckling = Inf(1,6);
            if F1 < 0, FOS_buckling(1) = Pcr_1 / abs(F1); end
            if F2 < 0, FOS_buckling(2) = Pcr_2 / abs(F2); end
            if F3 < 0, FOS_buckling(3) = Pcr_3 / abs(F3); end
            if F4 < 0, FOS_buckling(4) = Pcr_4 / abs(F4); end
            if F5 < 0, FOS_buckling(5) = Pcr_5 / abs(F5); end
            if F6 < 0, FOS_buckling(6) = Pcr_6 / abs(F6); end

            % Yield & Axial FOS
            % If sigma_comb_* is zero or negative, this will produce a warning or NaN; keep as-is
            FOS_yield = [YS/sigma_comb_1, YS/sigma_comb_2, YS/sigma_comb_3, ...
                         YS/sigma_comb_4, YS/sigma_comb_5, YS/sigma_comb_6];

            % Avoid division by zero for axial FOS: if any Fi == 0 make that fail (FOS zero)
            if any([F1,F2,F3,F4,F5,F6] == 0)
                FOS_axial = zeros(1,6);
            else
                FOS_axial = [abs(axial_max/F1), abs(axial_max/F2), abs(axial_max/F3), ...
                             abs(axial_max/F4), abs(axial_max/F5), abs(axial_max/F6)];
            end

            % Stability check: center-of-pressure must fall within plate radius for both triples
            s_margin = 0.02;
            r_plate = r_bot_act + space;
            reac_rad = sqrt(x_reac^2 + y_reac^2);
            pass_check = (reac_rad <= (r_plate - s_margin));

            % Criteria check: require all FOS >= 2 and pass_check true
            if any(FOS_axial < 2) || any(FOS_buckling < 2) || any(FOS_yield < 2) || ~pass_check
                anyYawFailed = true;
                break  % break yaw loop; this (phi,theta) fails full yaw test
            end

        end % for psi

        % After yaw loop: if ALL yaw passed, save this (phi,theta) as successful
        if ~anyYawFailed
            Results(idx).maxRoll = phi;
            Results(idx).MaxPitch = theta_DCM;
            Results(idx).FOS_yield = FOS_yield;
            Results(idx).FOS_axial = FOS_axial;
            Results(idx).FOS_buckling = FOS_buckling;
            Results(idx).F1 = F1;
            Results(idx).F2 = F2;
            Results(idx).F3 = F3;
            Results(idx).F4 = F4;
            Results(idx).F5 = F5;
            Results(idx).F6 = F6;
            Results(idx).Rot_act_global = Rot_act_global;
            Results(idx).Bus_CM_global = Bus_CM_global;

            % ---- Store actuator lengths and unit vectors ----
            Results(idx).L = [L1(:); L2(:)];                 % 6x1 vector of lengths
            Results(idx).Uvec = [u1 u2 u3 u4 u5 u6];        % 3x6 matrix of unit vectors

            idx = idx + 1;
        else
            % do nothing: this (phi,theta) did not survive a full 360 yaw test
            % continue to next pitch (we do NOT break here)
        end

    end % for theta_DCM
end % for phi

%======================= DISPLAY RESULTS ============================
all_FOS_yield    = vertcat(Results.FOS_yield);
all_FOS_buckling = vertcat(Results.FOS_buckling);
all_FOS_axial    = vertcat(Results.FOS_axial);

% Conservative (minimum) factors of safety across all configurations
minYieldFOS    = min(all_FOS_yield(:));
minAxialFOS    = min(all_FOS_axial(:));
minBucklingFOS = min(all_FOS_buckling(:));

% Store corresponding max safe orientation (last fully safe roll/pitch in Results)
maxRoll_deg  = max([Results.maxRoll]);
maxPitch_deg = max([Results.MaxPitch]);

% Summary table of conservative values
Summary = table(maxRoll_deg, maxPitch_deg, minYieldFOS, ...
                minAxialFOS, minBucklingFOS, omega, ...
                'VariableNames', {'Max Roll (deg)','Max Pitch (deg)', ...
                                  'Yield FOS','Axial FOS','Buckling FOS','omega (deg/s)'})

%======================= FORCE SUMMARY ==============================
all_F1 = [Results.F1];
all_F2 = [Results.F2];
all_F3 = [Results.F3];
all_F4 = [Results.F4];
all_F5 = [Results.F5];
all_F6 = [Results.F6];
all_forces = [all_F1(:); all_F2(:); all_F3(:); all_F4(:); all_F5(:); all_F6(:)];
[maxForce, maxIdx] = max(abs(all_forces)); % absolute max and its index

% Determine which configuration (and actuator) the max force came from
numResults = numel(Results);
actPerResult = 6;

% Find which result entry the max came from
resultIdx = ceil(maxIdx / actPerResult);
actuatorNum = mod(maxIdx-1, actPerResult) + 1; % 1..6

% Retrieve the corresponding roll, pitch, and CM of actuators
maxRoll_atForce  = Results(resultIdx).maxRoll;
maxPitch_atForce = Results(resultIdx).MaxPitch;
Satellite_CM = Results(resultIdx).Bus_CM_global;
Rot_CM = Results(resultIdx).Rot_act_global;

% ---- Retrieve actuator lengths and orientations for the worst-case result ----
L_vals = -2*U_joint + Results(resultIdx).L;        % 6x1 lengths
U = Results(resultIdx).Uvec;         % 3x6 unit vectors

% Compute actuator angles from vertical (degrees)
vertical = [0;0;1];
angles_deg = zeros(6,1);
for i = 1:6
    ui = U(:,i);
    angles_deg(i) = acosd( dot(ui, vertical) / (norm(ui)*norm(vertical)) );
end        

Output: