Defining Non Linear Vector Constraints in JuMP

Specific Domains Optimization (Mathematical)
1

I’m trying to implement pseudospectral methods using JuMP & Ipopt. The state is n-dimensional and control m-dimensional and the values are evaluated at N+1 collocation points; in this particular problem, n=3, m=3, and N = 120.

For pseudospectral methods, one constraint is that the numerical approximation of the derivative 2/TDX = f_all(X), where X is the stacked state [x1;x2…;x_{N+1}] and f_all = [f(x1); f(x2);…;f(x_{N+1})].

@NLconstraint(mod, f_all(dyn, state, N, n, m) - 2/Tp*D*state[1:n*(N+1)] == zeros(n*(N+1)))

dyn is the non-linear dynamics function and f_all returns a stacked vector of dyn() evaluated for each x.

dyn(p,w) = 0.25*((1-norm(p)^2)*w - 2*skew(w)*p + 2*dot(w,p)*p);
JuMP.register(mod, :dyn, 2, dyn, autodiff=true);

I keep running into numerous errors regardless of how I try to define the nonlinear functions:

ERROR: Cannot multiply a quadratic expression by a variable

I also tried breaking apart the n-D derivative into n different @NLexpression:

@NLexpression(mod, dyn1_expr[1:n], [0.25*(p[1]^2-p[2]^2-p[3]^2+1)*w[1] +  0.5*(p[1]*p[2]-p[3])*w[2] + 0.5*(p[2]-p[1]*p[3])*w[3]; 0.5*(p[1]*p[2]+p[3])*w[1] + 0.25*(p[2]^2 - p[1]^2-p[3]^2)*w[2] + 0.5*(p[2]*p[3]-p[1])*w[3]; 0.5*(p[1]*p[3]-p[2])*w[1] + 0.5*(p[1]+p[2]*p[3])*w[2] + 0.25*(p[3]^2-p[2]^2-p[1]^2+1)*w[3]])

and received:

ERROR: Unexpected vector JuMP.GenericQuadExpr{Float64,JuMP.Variable}[1.736231682755747e-5 state[1]² + 0.00033454250210412234 state[1]*state[2] - 1.736231682755747e-5 state[2]² - 0.0006879201991741926 state[1]*state[3] - 1.736231682755747e-5 state[3]² - 0.00033454250210412234 state[3] + 0.0006879201991741926 state[2] + 1.736231682755747e-5,-0.00016727125105206117 state[1]² + 3.472463365511494e-5 state[1]*state[2] + 0.00016727125105206117 state[2]² + 0.0006879201991741926 state[2]*state[3] - 0.00016727125105206117 state[3]² + 3.472463365511494e-5 state[3] - 0.0006879201991741926 state[1],-0.0003439600995870963 state[1]² - 0.0003439600995870963 state[2]² + 3.472463365511494e-5 state[1]*state[3] + 0.00033454250210412234 state[2]*state[3] + 0.0003439600995870963 state[3]² - 3.472463365511494e-5 state[2] + 0.00033454250210412234 state[1] + 0.0003439600995870963] in nonlinear expression. Nonlinear expressions may contain only scalar expressions.

Any idea how I could go about defining an n*(N+1) dimensional constraint by evaluating the dynamics f(x) N+1 times?

2 Likes
2

JuMP’s syntax and automatic differentiation utilities were not designed to handle this case. There’s a workaround for multidimensional input (Julia+JuMP: variable number of arguments to function - Stack Overflow) but not for multidimensional output. I’d recommend looking into tools like CasADi, or you can use ReverseDiff to compute derivatives and manually communicate them to Ipopt via MathProgBase.

3

If I were to break apart the dynamics function f into n separate functions, would that circumvent the issue? Since I know the matrix D from D*X, I could grab the i’th row and call a function f_i that also returns the i’th output of f(x). It seems inefficient, but that would allow me to pass in an n dimensional input and return one output at a time.

1 Like
4

Yes, you can do that. I won’t claim that it’s pretty though.

5

Is it still the case that matrix-vector products cannot be used in nonlinear constraints?

1 Like
6

@acauligi, If you don’t mind . If you got any working solution to this problem, kindly provide it here. Thank you.