Alright, getting a little bit closer, but still quite a few big problems. First,
@variable(m,T) creates a scalar variable called
T, and you’re currently trying to constrain each element of the matrix
T_ to be equal to
T. That’s probably not what you want.
Second, are you aware that JuMP’s variables are assumed to be real?
conj applied to a real matrix
A just returns
A. Are you looking for
lets go for a matrix exponential for now.
OK, in Julia, the
e.^A is used to signify a broadcasted operation (in this case, an elementwise-applied exponential function). If you were operating on a numeric array, you’d want
expm(A) for a matrix exponential (or just
exp(A) in the upcoming version of Julia).
Next, it’s important to note the difference between linear/quadratic modeling and nonlinear modeling. All of the expressions you construct up to
@constraint(m,Exp1.==e.^(*(T,L))) are linear or quadratic in the decision variables you create using
@variable. That code works because various operators (like
+) and functions are overloaded for e.g.
GenericAffExpr types. But JuMP’s nonlinear interface unfortunately doesn’t work quite the same way as its linear/quadratic interface. It’s quite a bit harder to use.
The following nonlinear syntax notes are especially relevant to your problem:
All expressions must be simple scalar operations.
This means that creating an expression that represents the matrix exponential of
T * L unfortunately can’t work, as it is not a scalar expression.
There is no operator overloading provided to build up nonlinear expressions
All nonlinear expressions must be inside of macros.
Before, you were trying to construct a nonlinear expression outside of a JuMP macro, which is not supported, so the code in your latest post is better in that sense, but still not quite right.
@expression is only for linear/quadratic modeling, so you’d need
@NLexpression, but again,
@NLexpression can only be used to construct scalar nonlinear expressions.
Here’s a slightly updated version of the first part of your problem formulation:
using Ipopt # just because that's what I have installed right now
C(A, B) = A * B - B * A
m = Model(solver=IpoptSolver())
N = 2
L = 100000
@variable(m, Rdiag[i=1:N]) # first create the variables that will fill the array
R = flipdim(Matrix(Diagonal(Rdiag)),2) # create the array with those variables along its diagonal
Q = Diagonal(Qdiag)
T = kron(Q, Diagonal(ones(2))) + kron(Diagonal(ones(2)), conj(Q)) + kron(R, conj(R))
LT = L * T
# since syntax notes state that
# > 'AffExpr and QuadExpr objects cannot currently be used inside nonlinear expressions'
# need to create matrix of slack variables `LT_` and set them equal to `LT`
@variable(m, LT_[1:size(LT, 1), 1:size(LT, 2)])
for i in eachindex(T)
@constraint(m, LT_[i] == LT[i])
The next part, involving the nonlinear expressions for
expm and for e.g.
kron(C(Q,R),C(conj(Q),conj(R))) (a matrix of quartic expressions) is the tricky bit, but unfortunately this is all I have time for right now. Some hints for you or others who want to help out:
- if you use the same slack variable trick as above for
C(Q, R) and
C(conj(Q),conj(R)), you can get
kron(C(Q,R),C(conj(Q),conj(R))) as a matrix of quadratic expressions, and then you can use the slack variable trick again to reduce this to a matrix of plain variables.
- since you can only use scalar nonlinear expressions, you’ll have to do the whole computation of the objective given
C(Q, R) and
C(conj(Q),conj(R)) in one function that you then register with JuMP. Nonlinear functions can’t take matrix or vector arguments though; search Discourse for some solutions to this problem.
- there’s currently no generic
expm function, so automatic differentiation won’t work; you’ll have to supply derivatives of the manually registered function from 2. yourself (not too hard in this case)