ControlSystems.jl icon indicating copy to clipboard operation
ControlSystems.jl copied to clipboard

Adjustment to dlqr and dlyap

Open ueliwechsler opened this issue 5 years ago • 2 comments

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?

ueliwechsler avatar Jun 11 '19 19:06 ueliwechsler

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.

baggepinnen avatar Jun 11 '19 19:06 baggepinnen

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.

mfalt avatar Jun 11 '19 20:06 mfalt