hi all, I’m very new to Julia. I’m following these lectures to solve optimization problem.
I will be providing a toy problem (2 DoF) which mimics my actual problem (12 DoF).
I have a set of nonlinear non-autonomous differential equations . (dx, dy). For some value of parameter “a” this system has a periodic solution, with a period “T”.
‘’'julia
function dzdt(x,y,t,a)
u = (1-(x^2)/4)-y^2
dx = -4*y+x*u + a*x*exp((10^-3)*t)
dy = x+y*u + a*y*exp((10^-2)*t)
[ dx, dy ]
end
‘’’
There might be many methods to solve these kind of problems. I’m using collocation and spectral methods for periodic system to find periodic solution of given system.
Key idea
-
let “x(t), y(t)” be periodic solution to system “dzdt” for some “a”. We can easy represent a some random periodic function in Fourier series. Instead of continuous time here we deal with discrete time.
1.a. If we can represent solutions-> ( x(t), y(t) ) in terms of some Fourier series → ( F1(t), F2(t) ).
1.b. We can take time derivative of these solutions → ( d (x(t), y(t) ) /dt ) which should be our given set of differential equations.
1.c. Difference between guess periodic solution (actually the derivative of the guess solution) and actual differential equations should be zero. Then we can say that we found the periodic solution. -
We construct this solution by assuming for some discrete time points (t = t1,t2,t3…tn) and at these time points ( x(t) , y(t) ) pass through some collocation points ( ( x1,y1), (x2,y2), (x3,y3),… (xn,yn) ). These (xi’s, yi’s) are our guess points.
-
We have matrix called "D Matrix for Fourier " . when this “D” matrix is pre-multiplied with the vector (v) it gives the approximate solution of the derivative of the (v) vector.
Here is the D Matrix
‘’'julia
function D(N::Integer)
h=2*pi/N;
ns=range(1,N-1,step=1);
col1(nps)=0.5*((-1)^nps)/sin(nps*h/2);
col=[0,col1.(ns)...];
row=[0,col[end:-1:2]...];
Toeplitz(col,row)
end
‘’’
Where “N” is the number of discrete time points.
Discrete points in time. Since I don’t know the period of system I will take some high value as final time “tfinal” . Here I will construct some initial guess solution points for x(t),y(t).
‘’’ julia
tfinal = 1.1*pi;
tpoints = collect( range(1,N,step=1) )*tfinal/N;
xguess = sin.( (2*pi/tfinal)*tpoints )*2.0;
yguess = -sin.( (2*pi/tfinal)*tpoints )*0.5;
‘’’
Here I’m writing a function to find the difference of assumed derivative solution and actual system.
I have written this function such that this can be used as a constraint but not as a objective function.
My objective function could be “0” . as of now.
‘’’ julia
function dxlist(xs,ys,tf,a)
ns = 2;
ts = collect( range(1,N,step=1) )*tf/N;
xytsZip = zip(xs,ys,ts);
dxD0 = [ dzdt(x,y,t,a) for (x,y,t) in xytsZip ];
dxD = reduce(hcat, dxD0)';
xlyl = reshape([xs;ys],N,ns);
dxF = (Dmat*xlyl)*(2.0*pi/tf);
err = dxD - dxF;
[ vcat(err'...) .-10^(-10) ; -vcat(err'...) .+10^(-10) ]
end
function cons(x)
# nps=(length(x[:])-1) ÷ 2;
tf=x[end-1];
a=x[end];
xs1=x[1:N];
ys1=x[N+1:2*N];
dxlist(xs1,ys1,tf,a)
end
‘’’
initial guess points as a vector
‘’’ julia
a0 = 10^-3;
x0 = vcat( [ xguess; yguess; [ tfinal , a0] ] );
xlower1=push!( -3.0*ones(2*N), pi );
xlower=push!(xlower1,-10^-3);
xupper1=push!(3.0*ones(2*N),1.5*pi);
xupper=push!(xupper,10^-3);
consLower=-ones(4*N)*Inf;
consUpper=zeros(4*N);
‘’’
Here im calling optimization functions
‘’'julia
obj(x)=0.0;
model=ADNLPModel(obj, x0, xlower, xupper, cons, consLower, consUpper)
output=ipopt(model)
xstar=output.solution
fstar=output.objective
‘’’
Packages used are
‘’'julia
using LinearAlgebra
using DiffEqOperators, ForwardDiff
using ApproxFun
using FFTW
using ToeplitzMatrices
using ModelingToolkit
using DifferentialEquations
using NLPModels
using ADNLPModels
using NLPModelsIpopt
using DCISolver
using JSOSolvers
‘’’
Whole code in one place.
‘’’ julia
using LinearAlgebra
using DiffEqOperators, ForwardDiff
using ApproxFun
using FFTW
using ToeplitzMatrices
using ModelingToolkit
using DifferentialEquations
using NLPModels
using ADNLPModels
using NLPModelsIpopt
using DCISolver
using JSOSolvers
# Number of collocation points
N=3
function Dmatrix(N::Integer)
h=2*pi/N;
ns=range(1,N-1,step=1);
col1(nps)=0.5*((-1)^nps)/sin(nps*h/2);
col=[0,col1.(ns)...];
row=[0,col[end:-1:2]...];
D=Toeplitz(col,row)
end
Dmat=Dmatrix(N);
function dzdt(x,y,t,a)
u=(1-(x^2)/4)-y^2;
dx=-4*y+x*u+a*x;
dy=x+y*u+a*y;
[dx,dy]
end
# initial guess
tfinal=1.1*pi;
tpoints=collect(range(1,N,step=1))*tfinal/N;
xguess=sin.((2*pi/tfinal)*tpoints)*2.0
yguess=-sin.((2*pi/tfinal)*tpoints)*0.5
function dxlist(xs,ys,tf,a)
nstates=2
ts=collect(range(1,N,step=1))*tf/N;
xytsZip=zip(xs,ys,ts);
dxD0=[dzdt(x,y,t,a) for (x,y,t) in xytsZip];
dxD=reduce(hcat, dxD0)';
xlyl=reshape([xs;ys],N,nstates);
dxF=(Dmat*xlyl)*(2.0*pi/tf);
err=dxD-dxF;
[vcat(err'...).-10^(-10);-vcat(err'...).+10^(-10)]
end
function cons(x)
tf=x[end-1];
a=x[end];
xs1=x[1:N];
ys1=x[N+1:2*N];
dxlist(xs1,ys1,tf,a)
end
a0=10^-3;
x0=vcat([xguess;yguess;[tfinal,a0]]);
obj(x)=0.0;
xlower1=push!(-3*ones(2*N),pi);
xlower=push!(xlower1,-10^-3)
xupper1=push!(3*ones(2*N),1.5*pi);
xupper=push!(xupper,10^-3)
consLower=-ones(4*N)*Inf;
consUpper=zeros(4*N)
# println("constraints vector = ",cons(x0))
model=ADNLPModel(obj,x0,xlower,xupper,cons,consLower,consUpper)
output=ipopt(model)
xstar=output.solution
fstar=output.objective
‘’’
When I implement this code i’m getting error. Any suggestion on how find orbit of such nonlinear non-autonomous system would be highly appreciated. Kindly help me in solving this toy problem. Thanking you all.