Julia求解带约束的非线性规划问题（改写Matlab程序）？

综合讨论区 性能优化
1

Matlab核心部分代码如下（可简化的不多了，索性都贴上了，核心是fmincon函数，需要读取grain_spot_list.txt文件https://1drv.ms/t/s!Alc9TKg71lAsmYACBjhpi6TyQe8tOQ?e=vJa4mx
直接运行test函数即可。由于类似操作要做几万遍，用Matlab实在太耗时，希望能用Julia实现并提高效率，奈何水平有限，特来此求助，多谢 :pray:

function [refined_ori, refined_pos, refined_error] = test
spot_list = importdata('grain_spot_list.txt');
B = [2.191936266240916  0   0
    0   2.191936266240916   0
    0                   0   2.191936266240916];
x0 = [0.717821405534626
    0.692917298758609
    0.067808899396516
    -0.447251137765496
    0.533571784789791
    -0.717821405534626
    -0.533571784789791
    0.484939962837988
    0.692917298758609
    0
    0
    0]';
ConstraintFunction = @constraintfun;
options = optimoptions('fmincon','Display','off');
LB = [-1 -1 -1 -1 -1 -1 -1 -1 -1 -0.5 -0.5 -0.5];
UB = [ 1  1  1  1  1  1  1  1  1  0.5  0.5  0.5];
[refined_grain,fval,~,~]  = ...
    fmincon(@(x) fitgrains(x,spot_list,B),x0,...
    [],[],[],[],LB,UB,ConstraintFunction,options);
refined_ori = reshape(refined_grain(1:9),3,3);
refined_pos = refined_grain(10:12);
refined_error = fval;
end

function angle_sum = fitgrains(x,spot_list,B)
U = reshape(x(1:9),3,3);
pos = x(10:12);
refined_Gvs = zeros(size(spot_list,1),3);
Gv_angle = zeros(size(spot_list,1),1);
for i = 1:size(spot_list,1)
    rot = spot_list(i,2);
    phi1 = rot*pi/180;
    PHI = 0;
    phi2 = 0;
    U11 = cos(phi1)*cos(phi2)-sin(phi1)*sin(phi2)*cos(PHI);
    U21 = sin(phi1)*cos(phi2)+cos(phi1)*sin(phi2)*cos(PHI);
    U31 = sin(phi2)*sin(PHI);
    U12 = -cos(phi1)*sin(phi2)-sin(phi1)*cos(phi2)*cos(PHI);
    U22 = -sin(phi1)*sin(phi2)+cos(phi1)*cos(phi2)*cos(PHI);
    U32 = cos(phi2)*sin(PHI);
    U13 = sin(phi1)*sin(PHI);
    U23 = -cos(phi1)*sin(PHI);
    U33 = cos(PHI);
    Omega = [U11 U12 U13; U21 U22 U23; U31 U32 U33];
    spotsWeightedCentroid = spot_list(i,6:7);
    pos_rot = Omega*pos';
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    dety0 = 512.4836;
    detz0 = 512.0230;
    pixelysize = 0.067176;
    pixelzsize = 0.067176;
    Rtilt = eye(3);
    Lsam2sou = [-100     0     0];
    Lsam2det = 100;
    K_in = [pos_rot(1) pos_rot(2) pos_rot(3)] - Lsam2sou;
    dety = spotsWeightedCentroid(1);
    detz = spotsWeightedCentroid(2);
    spot_pos = Rtilt*[Lsam2det (dety-dety0)*pixelysize (detz-detz0)*pixelzsize]';
    K_out = spot_pos - pos_rot;
    refined_Gvs(i,:) = normr(normr(K_out')-normr(K_in));
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Gv_angle(i) = acos(dot(normr(Omega*U*B*spot_list(i,8:10)'),normr(refined_Gvs(i,:))));
end
angle_sum = sum(abs(Gv_angle));
end

function [c, ceq]= constraintfun(x)
c = [];
U = reshape(x(1:9),3,3);
ceq = U*U'-eye(3);
end
2

确认是 matlab 慢，而不是你用的求解方法、求解器慢？

我理解的优化问题建模完毕后就看你选的求解器了。
不是自己写的求解器，你用 julia 也是调别人的求解器，不会太多的。

1 个赞
3

基本搞定了，速度比fmincon快一半左右，但只能用Ipopt实现（不能并行），用NLopt优化不动（thread-safe但直接返回初值）。我把程序简化了下重新开了个讨论，能不能帮忙看一下？