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/T*D*X = 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?