Solve cannot compute a matrix system

[DanielMartensson]

Hi!

Let's say that I want to compute X from the algebraic Riccati equation (ARE).

A' X + X A - X B inv(R) B' X + Q = 0

I just do this:

syms X

X = solve(A' X + X A - X B inv(R) B' X + Q == 0, X)

The results is this:

error: Python exception: AttributeError: MutableDenseMatrix has no attribute is_Relational. occurred at line 3 of the Python code block error: called from python_cmd at line 176 column 5 solve at line 119 column 9

A is 6x6, B is 6x1, R is a constant, Q is 6x6 identity matrix and X will going to be 6x6 matrix.

To solve an ARE equation, I use the care command from control package. But I want to have the skills to use the theory with symbolic tools, to compute the ARE equation. Not rest my knowledge on a built in function, which I don't know how it works, only use. https://octave.sourceforge.io/control/function/care.html

Information: GNU octave version: 4.0.3 Symbolic version: 2.4.0-2 Operative system: Lubuntu Linux 17.04

[cbm755]

Let's try to narrow this down a bit.

1. syms X makes symbolic variable which is usually (always?) thought of as a scalar.

2. Maybe you want something like X = sym('x', [6 6])?

3. Regarding the error message, can you come up with a simpler matrix problem which gives the same error message? I assume it is coming from == 0, but would be good to confirm.

[DanielMartensson]

Let's say I want to compute P from:

PA + A'P - PBR^-1B'P +Q = 0

This is the Riccati equation.

Where A, B, Q, R are know matrices. R sometimes is just a constant.

Example:

>> A
A =

   2   3
   2  -2

>> B
B =

   2   1
   0   1

>> Q
Q =

   5  -2
  -2   1

>> R
R =

   2   2
   2   4

>> X = sym('x', [2 2])
X = (sym 2×2 matrix)

  ⎡x₁₁  x₁₂⎤
  ⎢        ⎥
  ⎣x₂₁  x₂₂⎦

>> X = solve(A' * X + X * A - X * B * inv(R) * B' * X + Q == 0, X)
warning: Using rat() heuristics for double-precision input (is this what you wanted?)
warning: called from
    sym at line 225 column 7
    numeric_array_to_sym at line 36 column 14
    sym at line 203 column 7
    mtimes at line 65 column 5
warning: Using rat() heuristics for double-precision input (is this what you wanted?)
warning: Using rat() heuristics for double-precision input (is this what you wanted?)
warning: Using rat() heuristics for double-precision input (is this what you wanted?)
error: Python exception: AttributeError: MutableDenseMatrix has no attribute is_Relational.
    occurred at line 3 of the Python code block
error: called from
    python_cmd at line 176 column 5
    solve at line 119 column 9
>>

Not even this works:

>> X = solve(A' * X + X * A - X * B * inv(R) * B' * X + Q == [0;0], X)
warning: Using rat() heuristics for double-precision input (is this what you wanted?)
warning: called from
    sym at line 225 column 7
    numeric_array_to_sym at line 36 column 14
    sym at line 203 column 7
    mtimes at line 65 column 5
warning: Using rat() heuristics for double-precision input (is this what you wanted?)
warning: Using rat() heuristics for double-precision input (is this what you wanted?)
warning: Using rat() heuristics for double-precision input (is this what you wanted?)
error: Python exception: AssertionError
    occurred at line 8 of the Python code block
error: called from
    python_cmd at line 176 column 5
    binop_helper at line 60 column 5
    ineq_helper at line 35 column 5
    eq at line 91 column 5
error: evaluating argument list element number 1
>>

I have been trying to solve that with fsolve, but I won't work.

[cbm755]

I will try look at this later, but is this really the simplest problem you can find? Is it the transposes? Is it the inverse? It might be happening before u even call solve.

Please help narrow this down to a minimal example.

[mtmiller]

A much simpler example that results in the same error message:

x = sym ('x', [2, 1])
solve (x == 0)

But breaking the matrix of equations into individual arguments to solve works.

It looks like Matlab solve works this way, but ours does not yet support this. Example given on the Matlab public help page:

syms u v
eqns = [2*u + v == 0, u - v == 1];
S = solve(eqns, [u v])

also fails with the same error in Octave. However, solve(eqns(1), eqns(2), u, v) works.

[mtmiller]

From reading the Matlab docs, it doesn't look like solve(eqn1, eqn2, var1, var2) is supported at all. We could try to support both, but eliminating this form completely might help simplify the logic a lot.

Can someone with access to Matlab very whether this works?

syms x y
[X, Y] = solve(x*x == 4, x == 2*y, x, y)

Or does it have to be written like this?

[X, Y] = solve([x*x == 4, x == 2*y], [x, y])

[cbm755]

I'm reasonably sure SMT at least used to supported (I remember worrying about ambiguity so would not have implemented it this way out of my choice).

Worth checking if this is supported in Matlab SMT:

[X, Y] = solve({x*x == 4, x == 2*y}, {x, y})
[X, Y] = solve({x*x == 4, x == 2*y}, [x, y])

[mtmiller]

In Matlab 2017, the following forms are all equivalent:

syms x y
[X, Y] = solve (x*x == 4, x == 2*y, x, y)
X =

 -2
  2

Y =

 -1
  1

[X, Y] = solve (x*x == 4, x == 2*y, [x, y]);  %% same
[X, Y] = solve ([x*x == 4, x == 2*y], [x, y]);  %% same
[X, Y] = solve ([x*x == 4, x == 2*y], x, y);  %% same

@cbm755 the two forms you suggested with cell arrays are not supported.

[cbm755]

thanks, I will get to this eventually but I certainly don't mind if you want to do it!