Why is my step response going down instead of up?

I am trying to plot the step response of a 6th order ITAE system as part of a homework assignment. One of the criteria is to use a undamped natural frequency that is 2 times of one that i used in a previous part of the problem. When I plot the step response, the plot is going down, which seems off.

I am wondering if there is a problem with my code and am hoping somebody can point out my error. My code is segmented into part i), ii) and iii). Part iii) is where I deal with the 6th order ITAE.

Thanks in advance

% i) Calculate the undamped natural frequency and damping ratio required
% for a standard generic second-order system to achieve (assuming unit step input). 

% WHAT are the associated sytem eigenvalues? Plot the unit step response
% for the result and demonstrate how well the design criteria are
% met (normalize your output to ensure the final value is 1.0). 
% Display the resulting rise time, peak time, settling time, and percent
% overshoot on your graph.  USE Percent overshoot = 3, settling time=3

% Values used from chapter 1.
% Mass values 
m1 = 1; 
m2 = 2;
m3 = 3; 

% Spring Coefficients 
k1 = 1; 
k2 = 2; 
k3 = 3; 
k4 = 4; 

% Damping Coefficients 
c1 = 1; 
c2 = 2;
c3 = 3; 
c4 = 4;

% A, B, C and D matrix generated from chapter 1. 
A = [0 0 0 1 0 0; 
     0 0 0 0 1 0; 
     0 0 0 0 0 1; 
     ((-1)*((k1+k2)/m1)) (k2/m1) 0 ((-1)*((c1+c2)/m1)) (c2/m1) 0; 
     (k2/m2) ((-1)*((k2+k3)/m2)) (k3/m3) (c2/m2) ((-1)*((c2+c3)/m2)) (c3/m2); 
     0 (k3/m3) ((-1)*((k3+k4)/m3)) 0 (c3/m3) ((-1)*((c3+c4)/m3))];

B = [0 0 0;
     0 0 0; 
     0 0 0;
     (1/m1) 0 0;
     0 (1/m2) 0;
     0 0 (1/m3)];
 
 C = [1 0 0 0 0 0; 
      0 1 0 0 0 0; 
      0 0 1 0 0 0];
 
 D = [0 0 0; 
      0 0 0; 
      0 0 0];

sys_1 = ss(A,B,C,D);


% Specifying percent overshoot (3%) and settling time (3s): 
PO = 3; 
ts = 3; 

% Finding damping ratio and undamped natural frequency 
term = pi^2 + log(PO/100)^2; 
zeta = log(PO/100)/sqrt(term) % Damping ratio 
wn = 4/(zeta*ts) % Natural frequency from settling time and zeta.

% Setting up generic desired second order system 
num1 = wn^2;
den1 = [1 2*zeta*wn wn^2];

% Desired control law eignevalues 
DesEig1 = roots(den1) 
Des1 = tf(num1,den1) % Create desired systems from num2 and den2 

% Plotting 
figure();
td = [0:0.01:40.0];
step(Des1,td); 
% We will utilize this to display desired info. Commented out for now
stepinfo(Des1); 


%% Part ii) Augment desired second order system eigenvalues 
% Need 4 additional eigenvalues since this is a 6th order system 
% The first additional eigenvalue was chosen to be 10 times the
% real part of the eigenvalue obtained in part i) 

DesEig3 = real(DesEig1(1))*10; % We add the real() to isolate real part from complex number
DesEig4 = DesEig3-1;
DesEig5 = DesEig4-1;
DesEig6 = DesEig5-1;


%% Part iii) Computate the ITAE 6th-order coefficients & Eigenvalues.
% Undamped natural freq is 2 times that of part i). 
% Plot the 6th order ITAE and the augmented 6th order step response with 
% the dominant second order step response of part i). 

w_itae = wn*2
num_itae = (w_itae)^2; % NEEDS TO BE 2 TIMES w used in part i). Question is which w? or all? 


% 6th order ITAE 
den_itae = [1 3.25*wn 6.6*wn^2 8.6*wn^3 7.45*wn^4 3.95*wn^5 wn^6]; 

% Finding Eigenvalues from ITAE characteristic polynomial
DesEig_Itae = roots(den_itae) 
Itae_tf = tf(num_itae,den_itae);

figure();
td = [0:0.01:60.0];
step(Itae_tf,td);

Answer 1

[Pat Gipper]

I would focus on your equation for zeta on line 60. It is resulting in a negative value (negative damping=bad!) which gives unstable eigenvalues in the remainder of your code.

[Justin Goh]

Thank you so much. My textbook even had absolute values around the zeta.