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

Matlab核心部分代码如下（可简化的不多了，索性都贴上了，核心是fmincon函数，需要读取grain_spot_list.txt文件https://1drv.ms/t/s!Alc9TKg71lAsmYACBjhpi6TyQe8tOQ?e=vJa4mx

``````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

``````

1 个赞