ControlSystems.jl
ControlSystems.jl copied to clipboard
Adjustment to dlqr and dlyap
I came across some peculiar behaviour using dlyap
and dlqr
.
To my knowledge, the result from dlyap
should always return a symmetric matrix, however, when running the code below. The Matrix P
does not pass the issymmetric
test. Additionally, I would suggest, using broadcasting in the definition of dlqr
The first part is to construct a stable system dynamcis As
using LinearAlgebra
using ControlSystems
A = [1 -2; 0 1]
B = Matrix([1 1]')
R = 1
Q = [2 0; 0 2]
F = dlqr(A,B,Q,R)
This returns an error since there is no broadcasting in the function definition
(I would suggest using broadcasting in dlqr
K = (B'*S*B .+ R)\(B'S*A)
.
We can fix this, by using the following.
R = ones(1,1)
F = dlqr(A,B,Q,R)
As = A - B*F
P = dlyap(As, Q)
issymmetric(P)
Now, I want to use the resulting P
, but since it is not symmetric, as an example, I cannot use it for the construction of an Ellipsoid directly, I have to call Symmetric(P)
first.
Is this behaviour wanted or is this a bug?
If we use . +R
, we can not allow the user to send in uniform scaling I
, as that broadcasts as the matrix of ones. Another method should paobably be added for scalar, or use the internal to_mat
The non symmetric result is not intended behavior, its probably numerical noise.
I agree, we really want the behavior of the operation +
, it is bad if we accidentally accept vectors where we expect 2x2 for example. However, in the case of a scalar, we could probably handle it as baggepinnen says.
I tried the symmetry on my computer, and you are right that it is not isymmetric
, but (P-P')[2]≈5e-15
, so it seems to just be numerical noise. One option is to output Symmetric(P) from the funktion, which basically masks the lower half, but I'm not sure if this would crate mpre problems.
Also note that there is no need to make B
into a Matrix in your first case, the vector of right shape B = [0,1]
should work fine.