How to use Manifold projection in differential equations jl

Specific Domains
, ,
1

Hi all. I am trying to simulate a geometric control which we proposed in our paper. Since it involves a differential equation evolving on a manifold. I thought of using manifold projection as mentioned in Callback Library · DifferentialEquations.jl . However, I am getting instability in solving the equation. Also not able to set up the NLSOLVEJL_SETUP(). I am including my code here. It would be really helpful if anyone could suggest what I could do differently here.

----------------------------------CODE--------------------------------------

using DifferentialEquations;
using OrdinaryDiffEq;
using Sundials;
using LinearAlgebra;
using StaticArrays;
using ODEInterface;
using ODEInterfaceDiffEq;
using LSODA;

function toFindLogRot!(RotMat)

trace=tr(RotMat);
println(trace)
arg = (trace-1)/2;
if arg >= 1
    arg = 1
end
theta=acos(arg);
if  theta != 0
    omega=(1/(2*sin(theta)))*[RotMat[3,2]-RotMat[2,3];RotMat[1,3]-RotMat[3,1];RotMat[2,1]-RotMat[1,2]];
else
    omega=[1;0;0];
end 
return [omega' theta]'

end

----------------Manifold Part---------------
Since the rotation matrix has to satisfy transpose(R)*R =I, Can I write Like this??
function g(resid,R,p,t)

R = reshape(R,(3,3))
dR = transpose(R)*R - I
resid = reshape(dR,(1,:))

end

cb = ManifoldProjection(g,nlsolve=NLSOLVEJL_SETUP()) # How to setup nlsolve setup ??

function geometricControl!(dR,R,p,t)

R_if = reshape(R,(3,3))
p_if =[100*sin(t/20),100*sin(t/40),10];
p_if_dot = [5*cos(t/20);2.5*cos(t/40);0]
p_if_dotdot =  [-0.25*sin(t/20);-(25/400)*sin(t/40);0]
p_fr = transpose(R_if)*(p_ir-p_if);
r = sqrt(p_fr[1]^2+p_fr[2]^2+p_fr[3]^2)
theta=acos( p_fr[3] / r);
phi=atan(p_fr[2] , p_fr[1]);
Rd = [ cos(phi)*cos(theta) cos(theta)*sin(phi) -sin(theta);
            -sin(phi)            cos(phi)           0;
            cos(phi)*sin(theta) sin(phi)*sin(theta)  cos(theta)];
R_if_dot = reshape(dR,(3,3))
p_fr_dot = transpose(R_if_dot)*(p_ir-p_if)-transpose(R_if)*p_if_dot 
theta_dot = (-1/sqrt(p_fr[1]^2+p_fr[2]^2))*((r^2*p_fr_dot[3]-p_fr[3]*(p_fr[1]+p_fr[2]+p_fr[3]))/r^2)
phi_dot = (p_fr[1]*p_fr_dot[2]-p_fr[2]*p_fr_dot[1])/(sqrt(p_fr[1]^2+p_fr[2]^2))
Rd_dot = [-sin(phi)*cos(theta)*phi_dot-cos(phi)*sin(theta)*theta_dot cos(theta)*cos(phi)*phi_dot- 
           sin(theta)*sin(phi)*theta_dot -cos(theta)*theta_dot;
            -cos(phi)*phi_dot -sin(phi)*phi_dot 0;
            -sin(theta)*sin(phi)*phi_dot+cos(theta)*cos(phi)*theta_dot 
           sin(theta)*cos(phi)*phi_dot+cos(theta)*sin(phi)*theta_dot -sin(theta)*theta_dot]
er = toFindLogRot!(R_if*transpose(Rd))
eR = [0 -er[3] er[2];er[3] 0 -er[1];-er[2] er[1] 0];
omega = -beta * transpose(Rd)*er[4]*eR*Rd - transpose(Rd_dot)*Rd
R_if_dot = R_if*omega
dR .= reshape(R_if_dot,(1,:))

end

p=0
beta = 0.05
p_ir = [20,10,100];
tfinal = 2;
time_span = (0.0,tfinal);
R_if_0 = reshape(Matrix{Float64}(I,3,3),(1,:))
prob_di= ODEProblem(geometricControl!,R_if_0,time_span,p);
@time sol_di = solve(prob_di,CVODE_BDF(),callback=cb);

---------------------------ERROR----------------------

  1. Instability detected if nsolve is not used
  2. NLSOLVEJL_SETUP not defined

What I have tried —>
I have gone through this website Specifying (Non)Linear Solvers · DifferentialEquations.jl.
But code through an error CS and AD not defined (CS being chunk size, AD is Autodifferentiation (not sure about what AD mean)).

If possible please suggest what could be done about this.

2

Don’t pass in the default arg if you don’t need to. Just do cb = ManifoldProjection(g). See the example: https://diffeq.sciml.ai/stable/features/callback_library/#Example

3

Hey, I checked that out. But if I use manifold projection without NLSOLVEJL_SETUP I get “instability detected error”. I want to know how do I set up this NLSOLVEJL_SETUP. I could not understand how exactly it was created in callback documents on NLSOLVEJL_SETUP.

4

Instability detected likely is not due to the nonlinear solver that is chosen. I would check why the ODE is diverging itself.

But if you want to see how the nlsolver is defined, you can look at an example like:

This is something set to change (simplify) soon, and instead just take a NonlinearSolve option.