Julia中非线性向量约束的定义
提问于 2022-04-18 04:56:41
我试图最小化一个函数，它以一个向量作为输入，并且受到一些非线性约束。我对茱莉亚很陌生。我正在尝试用Ipopt.My实现伪谱方法，它是最优的，我使用的是代价函数和约束的梯度。像"ForwardDiff，ReverseDiff“这样的函数无助于找到向量函数的梯度。我发现@acauligi也面临着类似的问题。到目前为止我还没有找到任何解决办法。
using LinearAlgebra, DiffEqOperators, ForwardDiff, ApproxFun, FFTW, ToeplitzMatrices
using ModelingToolkit,DifferentialEquations,NLPModels,ADNLPModels,NLPModelsIpopt
using DCISolver,JSOSolvers
# Number of collocation points
N=31 # This number can go up to 200
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; backend =
ADNLPModels.ReverseDiffAD)
output=ipopt(model)
xstar=output.solution
fstar=output.objective在MatLab中，我在3分钟内得到了这个问题的解决方案(这个问题的解决方案是。当a=0时，系统的时间周期是π。)
我希望我能在朱莉娅身上更快地得到同样的结果。到目前为止，我在朱莉娅话语上已经问过了，我有任何建议。任何关于如何解决这个问题的建议都受到高度赞赏。谢谢你们所有人。
回答 1
Stack Overflow用户
回答已采纳
发布于 2022-09-28 11:57:18
我觉得你的代码有两个问题。第一,
xupper1=push!(3*ones(2*N),1.5*pi);
xupper=push!(xupper1,10^-3)然后，由于某种原因，另一个矩阵的Toeplitz矩阵的乘积与自动微分产生了误差。然而，以下工作：
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=vcat((Dmat*xlyl[:,1])*(2.0*pi/tf), (Dmat*xlyl[:,2])*(2.0*pi/tf));
err=vcat(dxD...) - dxF;
[err.-10^(-10);-err.+10^(-10)]
end最后，Ipopt返回正确的结果。
model=ADNLPModel(obj,x0,xlower,xupper,cons,consLower,consUpper)
output=ipopt(model)
xstar=output.solution
fstar=output.objective我还注意到使用Percival.jl更快
using Percival
output=percival(model, max_time = 300.0)
xstar=output.solution
fstar=output.objective请注意，ADNLPModels.jl正受到一些关注，并将显著改善。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/71912111
复制相关文章
腾讯云开发者
