I have the following code in Matlab that creates a stable system using a simple PID controller.
%% Define system
sys = ss(A,B,C,D);
%% New state space
sys_new = sys([1379],13);
%% Define s domain
Kp = 0;
Ki = 0; % No integrator
Kd = 278.7369;
Tf = 0;
sys_C = pid(Kp,Ki,Kd,Tf);
%% Controller design
sys_cl = feedback(sys_new,sys_C);
step(sys_cl)
which produces the following step response
The Julia version (or at least as close as i can make it is
using ControlSystems
sys = ss(A,B,C,D)
sys_new = sminreal(sys[1379,13])
Kp = 0
Ki = 0
Kd = 278.7369
Tf = 0
sys_C = pid(Kp,Ki,Kd; Tf=0.001, form=:parallel, state_space=true)
# notice the Tf has to be set to be non-zero here but not in the matlab version
sys_cl = feedback(sys_new,sys_C);
y,t = step(sys_cl);
plot(t,y[1,:])
This produces nonsense
The main difference I see in the Matlab version and the Julia version is that in the Julia version is that the PID controller needs to be converted into state_space and that requires the Tf be set to something non-zero which causes the step response to go bad. Indeed setting Tf=0.001 in the matlab version produces a similar looking response.
My question is how can I make the julia and matlab version agree?
Even with all computations performed using big floats, I still get that the unstable pole in the plant remains unstable after closing the loop:
using ControlSystems, BSON, LinearAlgebra
function LinearAlgebra.eigvals!(R::Matrix{BigFloat}, A::Matrix{BigFloat}; sortby::Nothing)
eigvals(A)
end
data = BSON.load("/tmp/state_space.bson")
A,B,C,D = getindex.(Ref(data), [:A, :B, :C, :D])
sys = ss(big.(A),big.(B),big.(C),big.(D))
# sys = ss(A,B,C,D)
sys_new = sminreal(sys[1379,13])
Kp = 0
Ki = 0
Kd = big(278.7369)
Tf = 0
sys_C = pid(Kp,Ki,Kd; Tf=1e-16, form=:parallel, state_space=true)
# notice the Tf has to be set to be non-zero here but not in the matlab version
L = sys_new*sys_C
sys_cl = feedback(sys_new, sys_C)
isstable(sys_cl)
I gave the matlab code in my first post and the step response flattens out unlike the julia version. What I don’t understand is how matlab is adding the feedback to the state-space model if the pid controller is in s-function form and not in state-space form. Is matlab taking the state-space system and turning it into a transfer function form?
I can’t run matlab, you can check what type sys_cl has, I’d expect them to convert to a descriptor system.
You can directly compute the poles of the closed loop system, that will tell you if you have any unstable ones left, before computing the step response.
The output type for sys_cl is a state_space system. Interestingly if I try to recreate the state_space model from the sys_cl ABCD I get an unstable system. For example
Are you doing the same thing in MATLAB and Julia? Above is what you do in MATLAB (unless you have a typo) – I don’t see any minimal realization there. Contrast that to what you do in Julia:
sys_new = sminreal(sys[1379,13])
So – in Julia, do you pick out row 1379 and column 13 of a (huge) input-output system, and do sminreal on that system? What are you doing in MATLAB?
[Sorry – my MATLAB CST competence is somewhat rusty…]
Yeah It does have a E matrix but its just the Identity matrix so I wouldn’t expect it to change anything. That being said if I do
%% matlab
sys_cp = ss(sys_cl.A, sys_cl.B, sys_cl.C,sys_cl.D)
sys_cp.E = sys_cl.E
step(sys_cl) # this gives the expected response.
According to the matlab docs E is
E — Matrix for implicit state-space models [] (default) | Nx-by-Nx matrix
Matrix for implicit or descriptor state-space models, specified as a Nx-by-Nx matrix. E is empty by default, meaning that the state equation is explicit. To specify an implicit state equation Edx/dt = Ax + Bu, set this property to a square matrix of the same size as A. See dss for more information about creating descriptor state-space models.
DescriptorSystems.jl allows system with non-identity E matrices. Still, it seems strange that the statespace system where the pid controller is filtered with an extremely large bandwidth wouldn’t yield the same result. I wonder if mstlab performs some model reduction in any of the steps.
No the original system without feedback control does not use E. Somehow adding PID feedback causes the combined system to have a non-identity E matrix i.e. a Descriptor System. How this happens? I’m not sure but the change does appear to make the system more numerically stable.
I am still puzzled by the very model of the system. The full model sys is of order 163 (the matrix A is 163x163). What is then the point in doing this?
Both in Matlab and in Julia I get (unsurprisingly) an error.
A pure derivative term, s, is not a proper system, it cannot be represented as a statespace system. This is why you have to add the low-pass filter in the pid constructor in order to be able to convert it to statespace. Non-proper systems can, however, be represented as descriptor systems:
I performed all computations involved in your problem using DescriptorSystems.jl and still got an unstable system, so if I were you I would double check the matlab results one more time.
