Cannot find a suitable controller for my state space system

%Define the parameters and system matrices as before
global g c k L m M Io f tau wn Ig K wb q z alpha ct

set(0,'defaulttextinterpreter','latex')
g=1.62; %lunar gravity
tau=-0.005; %envelope function 1/time constant
wn=117; %natural frequency of regolith foundation
m=2187500; %regolith foundation mass
M=3000;   %telescope mass 
Ig=1500;  %Moment of Inertia for telescope Used Ig=1/2*M*R^2 assumed radius of 1m
Io=(M*L^2)+Ig;
c=42485.625; %damping coefficient 
L=10;         %length of telescope 
k=29.96*(10^9); %spring constant of regolith
K=7*(10^10);
wb=18.85; %angular velocity of seismic input
q=M+m;
z=(M^2)*(L^2);
alpha=Io*(q)-z;
ct=2*sqrt(M*K)*0.1; %damping of telescope
% Define the state-space matrices
A = [0, 1, 0, 0; (-Io * k / alpha), (-Io * c / alpha), ((M * L * K - z * g) / alpha), M*L*ct/alpha; 0, 0, 0, 1; ((M * L * k) / alpha), ((M * L * c) / alpha), (((M * L * g *q) - (K*q)) / alpha), -q*ct/alpha];
B = [0; (Io / alpha); 0; (-M * L / alpha)];
Cgood=[1 0 0 0; 0 0 1 0]; %want to measure x and theta
C = [1, 0, 0, 0]; %x
C2= [0, 1, 0, 0];  %xdot
C3= [0, 0, 1, 0];  %theta
C4= [0, 0, 0, 1]; %thetadot
D = [0];
% Define the time vector
t = linspace(0, 1000, 1000000);
% define seismic input u 
u = @(time) (k * exp(tau * time) * sin(wb * time) + c * (tau * exp(tau * time) * sin(wb * time) + wb * exp(tau * time) * cos(wb * time)));
% Evaluate the input function for the entire time vector (so i can use u in lsim function)
u_matrix = arrayfun(u, t);

I am trying to design a controller for my state space system which is modelling a response of an inverted pendulum (telescope) on a base that is subject to ground motion (y). When making the state space model, I have u=ky+cydot where I define y as a damped sine wave to model an earthquake. However, whether it is using pid tuner, pole placement, or lqr control, I cannot design a controller that reduces oscillations of the system. The system is controllable and observable as the ranks of the controllability and observability matrices are 4. However, no matter what kind of technique I use, the controllers I design often do not change the sytem. I was wondering if there is something inherently wrong in the way I am approaching this problem. Thank you

Answer 1

[Sam Chak]

Hi @Ethan Nobre

The performance requirements are not clearly specified, other than addressing unwanted oscillations. Furthermore, it remains unclear how you utilized the PID tuner, designed the pole placement, or optimized the LQR controller. Generally, these tools prove to be less useful if appropriate requirements are not provided, unless one opts for the "I'm Feeling Lucky" approach. If the optimization-based LQR method is preferred, the cost function must be cleverly customized, as the conventional quadratic function may not effectively mitigate intense oscillations.

In the presented example, I devised a state-feedback controller with integral action using the pole placement method. This approach successfully eliminated unwanted oscillations, which I think is crucial in space applications as oscillations could potentially excite unmodeled vibration modes within the moonquake loading system. The configuration closely resembles the discrete-time control block diagram available at the following link:

https://www.mathworks.com/help/sps/ref/statefeedbackcontroller.html

%% Parameters
g       =  1.62;                % lunar gravity
tau     = -0.005;               % envelope function 1/time constant
wn      =  117;                 % natural frequency of regolith foundation
m       =  2187500;             % regolith foundation mass
M       =  3000;                % telescope mass 
Ig      =  1500;                % Moment of Inertia for telescope Used Ig=1/2*M*R^2 assumed radius of 1m
L       =  10;                  % length of telescope 
Io      = (M*L^2) + Ig;
c       =  42485.625;           % damping coefficient 
k       =  29.96e9;             % spring constant of regolith
K       =  7e10;        
wb      =  18.85;               % angular velocity of seismic input
q       =  M + m;
z       = (M^2)*(L^2);
alpha   =  Io*q - z;
ct      =  2*sqrt(M*K)*0.1;     % damping of telescope

%% Moonquake Loading System (4th-order)
A       =  [0, 1, 0, 0; -Io*k/alpha, -Io*c/alpha, (M*L*K - z*g)/alpha, M*L*ct/alpha; 0, 0, 0, 1; M*L*k/alpha, M*L*c/alpha, (M*L*g*q - K*q)/alpha, -q*ct/alpha];
B       =  [0; Io/alpha; 0; -M*L/alpha];
C       =  [1 0 0 0];           % x1 = x, x2 = xdot, x3 = theta, x4 = thetadot
D       =   0*C*B;
sys     = ss(A, B, C, D);
figure(1)
step(sys), grid on              % magnitude of uncompensated x is infinitesimally small

%% Pole-placement design (or use Ackermann's method)
Aa      =  [A zeros(4,1); -C 0] % extended state matrix (with x5 is integral of tracking error)

Aa = 5×5

1.0e+05 * 0 0.0000 0 0 0 -0.1370 -0.0000 0.0318 0.0000 0 0 0 0 0.0000 0 0.0136 0.0000 -2.3249 -0.0001 0 -0.0000 0 0 0 0

Bb      =  [B; 0];              % extended input matrix
wn      =   2.2136;             % natural frequency <-- only tune this parameter
poly5   =  [1 2.7*wn 4.9*wn^2 5.4*wn^3 3.45*wn^4 wn^5]; % don't alter unless you are expert
p       =  roots(poly5)         % desired poles

p = 
  -0.7647 + 2.9075i
  -0.7647 - 2.9075i
  -1.3939 + 1.2651i
  -1.3939 - 1.2651i
  -1.6595 + 0.0000i

Kk      =  place(Aa, Bb, p);    % use acker() if there are repeated poles

%% State-feedback controller with integral action
K       =  Kk(1:4)              % Full state-feedback gain matrix

K = 1×4

1.0e+12 * -0.0300 -0.0000 5.1106 0.0001

Ki      = -Kk(end)              % Tracking error integral gain

Ki = 500.7684

%% Compensated system
Acom    =  [A-B*K B*Ki; -C 0];
Bcom    =  [0; 0; 0; 0; 1];
Ccom    =  [C 0];
Dcom    =   D*Ccom*Bcom;
csys    = ss(Acom, Bcom, Ccom, Dcom);
S       = stepinfo(csys)        % check if step response characteristics are satisfied

S = struct with fields:
         RiseTime: 1.2682
    TransientTime: 2.5002
     SettlingTime: 2.5002
      SettlingMin: 0.9095
      SettlingMax: 1.0101
        Overshoot: 1.0096
       Undershoot: 0
             Peak: 1.0101
         PeakTime: 4.4955

%% Step responses
t       = 0:0.01:10;
% [y,x,t] = step(Acom, Bcom, Ccom, Dcom, 1, t);
[y,t,x] = step(csys, t);
x1      = [1 0 0 0 0]*x';   % x
x2      = [0 1 0 0 0]*x';   % x'
x3      = [0 0 1 0 0]*x';   % θ
x4      = [0 0 0 1 0]*x';   % θ'

%% Plot results
figure(2)
subplot(2,2,1)
plot(t, x1), grid on
xlabel('t'), ylabel('x_{1}'), title('x', 'fontsize', 14)

subplot(2,2,2)
plot(t, x2), grid on
xlabel('t'), ylabel('x_{2}'), title('x^{\prime}', 'fontsize', 14)

subplot(2,2,3)
plot(t, x3), grid on
xlabel('t'), ylabel('x_{3}'), title('\theta', 'fontsize', 14)

subplot(2,2,4)
plot(t, x4), grid on
xlabel('t'), ylabel('x_{4}'), title('\theta^{\prime}', 'fontsize', 14)

[Ethan Nobre]

Does this control system also work with the input being u where u is defined as

u = @(time) (k * exp(tau * time) * sin(wb * time) + c * (tau * exp(tau * time) * sin(wb * time) + wb * exp(tau * time) * cos(wb * time)));

[Sam Chak]

It didn't check this because I don't have full information about what this signal is or how you designed this time-dependent signal. Are you sure that you want to inject this highly oscillatory signal into the system?

This signal also appears to resemble some kind of seismic activity. Is this considered an external disturbance? If so, then you need to edit the input matrix B.

u    = @(time) (k*exp(tau*time).*sin(wb*time) + c*(tau*exp(tau*time).*sin(wb*time) + wb*exp(tau*time).*cos(wb*time)));
time = 0:0.005:1500;
plot(time, u(time)), grid on

[Ethan Nobre]

Yes this signal is supposed to model seismic activity and it is an external disturbance. I do not know how this would affect the state space matrices or the control system. When using lsim with the seismic signal as the input I have noticed that the gains become very high when using a control loop vs. when i do not have a controller, the system oscillates more frequently, but at a much lower amplitude.

tnew = linspace(0, 1000, 1000000);
u = @(time) (k * exp(tau * time) * sin(wb * time) + c * (tau * exp(tau * time) * sin(wb * time) + wb * exp(tau * time) * cos(wb * time)));
u_matrix = arrayfun(u, tnew);
[ynew,tnew,xnew] = lsim(csys,u_matrix', tnew);
x1new      = [1 0 0 0 0]*xnew';   % x
x2new      = [0 1 0 0 0]*xnew';   % x'
x3new      = [0 0 1 0 0]*xnew';   % θ
x4new      = [0 0 0 1 0]*xnew';   % θ'
figure
subplot(2,2,1)
plot(tnew, x1new), grid on
xlabel('t'), ylabel('x_{1}'), title('x', 'fontsize', 14)
xlim([0 10])
subplot(2,2,2)
plot(tnew, x2new), grid on
xlabel('t'), ylabel('x_{2}'), title('x^{\prime}', 'fontsize', 14)
xlim([0 10])
subplot(2,2,3)
plot(tnew, x3new), grid on
xlabel('t'), ylabel('x_{3}'), title('\theta', 'fontsize', 14)
xlim([0 10])
subplot(2,2,4)
plot(tnew, x4new), grid on
xlabel('t'), ylabel('x_{4}'), title('\theta^{\prime}', 'fontsize', 14)
xlim([0 10])

[Ethan Nobre]

@Sam Chak Also, how did you get that polynomial and why is wn tunable

[Sam Chak]

Hi @Ethan Nobre

I'm not a lunar expert, but I observed that the ordinary differential equations (ODEs) in the image are incomplete, as they do not indicate the channel through which the external disturbance enters. In your provided simulation, even though the complete code is not shown, I suspect that you introduced the disturbance into the same channel as the unit step function (desired input channel) in my simulation for the Compensated system. Is that indeed the case?

I recommend sketching the control block diagram and conducting simulations in Simulink. This approach will likely provide a more in-depth understanding of the Moonquake Loading System. The following simulation depicts the responses of x and theta in the Uncompensated system, influenced by the disturbing seismic force.

%% Parameters
g       =  1.62;                % lunar gravity
tau     = -0.005;               % envelope function 1/time constant
wn      =  117;                 % natural frequency of regolith foundation
m       =  2187500;             % regolith foundation mass
M       =  3000;                % telescope mass 
Ig      =  1500;                % Moment of Inertia for telescope Used Ig=1/2*M*R^2 assumed radius of 1m
L       =  10;                  % length of telescope 
Io      = (M*L^2) + Ig;
c       =  42485.625;           % damping coefficient 
k       =  29.96e9;             % spring constant of regolith
K       =  7e10;        
wb      =  18.85;               % angular velocity of seismic input
q       =  M + m;
z       = (M^2)*(L^2);
alpha   =  Io*q - z;
ct      =  2*sqrt(M*K)*0.1;     % damping of telescope

%% Moonquake Loading System (4th-order)
A       =  [0, 1, 0, 0; -Io*k/alpha, -Io*c/alpha, (M*L*K - z*g)/alpha, M*L*ct/alpha; 0, 0, 0, 1; M*L*k/alpha, M*L*c/alpha, (M*L*g*q - K*q)/alpha, -q*ct/alpha];
B       =  [0; Io/alpha; 0; -M*L/alpha];
C       =  [1 0 0 0];           % x1 = x, x2 = xdot, x3 = theta, x4 = thetadot
D       =   0*C*B;
sys     = ss(A, B, C, D);
ToF     = isstable(sys)         % True (1) or False (0)

ToF = logical
   1

% Responses to Seismic force
t       = 0:1e-3:1000;
u_seis  = (k*exp(tau*t).*sin(wb*t) + c*(tau*exp(tau*t).*sin(wb*t) + wb*exp(tau*t).*cos(wb*t)));
[y,t,x] = lsim(sys, u_seis, t); % x
theta   = [0 0 1 0]*x';         % θ

subplot(311)
plot(t, u_seis), grid on
title('Seismic force, u_{seis}'), ylabel('Force / Nm')
subplot(312)
plot(t, y), grid on
title('Positional response, x'), ylabel('Displacement / m')
subplot(313)
plot(t, theta), grid on
title('Angular response, \theta'), ylabel('Angle / rad'), xlabel('t / sec')

[Sam Chak]

Hi @Ethan Nobre

Also, how did you get that polynomial and why is wn tunable.

Technically, the pole placement method should be referred to as the Hurwitz Polynomial Placement method. The desired poles, as encountered in exercises found in most textbooks, do not magically manifest on their own. Rather, these poles are derived from solving the stable characteristic polynomial, whose roots reside in the left half-plane of the complex plane. In essence, control practitioners must acquire the skill of designing the Hurwitz polynomial.

However, this topic is not typically addressed in most undergraduate textbooks because it entails employing optimization strategies to determine the coefficients of the stable polynomial that showcases the desired time response characteristics. This may be covered in a postgraduate course such as "Optimal Control," usually offered by the Department of Mathematics (I have not taken this course yet).

Some books do provide the coefficients of the Hurwitz polynomials that satisfy certain integral performance indices, such as ISE, IAE, ITSE, ITAE. However, I personally prefer the Deadbeat performance criteria. The 5th-order Hurwitz Polynomial in my design can be obtained from "Modern Control Systems" by Dorf & Bishop.

As for why is tunable, you can observe that there is only one parameter in the Hurwitz Polynomial equation. One simply needs to manually tune a single parameter instead of relying on Genetic Algorithm (GA), Particle Swarm Optimization (PSO), or Reinforcement Learning (RL) to blindly search for the relatively high-dimensional set of 5 stabilizing control gains (in the Moonquake system).