目录

一、第1问

1.1 第一次遍历寻找浮标吃水深度MATLAB程序

1.2 第二次遍历寻找浮标吃水深度MATLAB程序

1.3 代入遍历结果浮标吃水深度MATLAB程序

二、第2问

2.1 遍历寻找满足条件重物球质量最小解MATLAB程序

2.2 加入多目标优化遍历寻找满足条件重物球质量MATLAB程序

2.3 代入遍历优化结果重物球质量MATLAB程序

三、第3问

3.1 遍历不同锚链型号与长度以及重物球的质量MATLAB程序

3.2 三次遍历代入锚链型号与长度以及重物球的质量求浮标吃水深度MATLAB程序

3.2.1 第一次遍历

3.2.2 第二次遍历

3.2.3 第三次遍历

3.3 代入遍历所得浮标吃水深度以及锚链型号与长度以及重物球的质量MATLAB程序

一、第1问

针对问题一，选取系泊系统相对静止状态为研究对象，进行静力学分析建立系泊系统各部分的受力平衡和力矩平衡模型，将所处海域的水深大小对系泊系统的垂直高度起限制作用作为系统纵向约束条件。将浮标的游动区域转化为游动半径进行研究，游动半径为浮标与锚的水平距离。假定浮标的吃水深度，对系泊系统的受力情况进行递推，得到系统各部分倾斜角度。对浮标的吃水深度在 0~2 米的范围变步长划分，进行循环遍历，建立垂直高度模型，寻找其满足系统纵向约束条件的值，并得到在不同风速下系泊系统的状态。结果发现，风速为12m/s 时，钢桶的角度为1.200 钢管的角度依次为1.158 、1.166 、1.174 、1.182，浮标的吃水深度为 0.683 米时，浮标的游动半径为 14.604 米，游动区域的面积大致为 670.029 平方米。风速为 24m/s 时，钢桶的角度为 4.563，钢管的角度依次为 4.410 、 4.438 、4.467 、4.496，浮标的吃水深度为 0.697 米，浮标的游动半径为17.780 米，游动区域的面积大致为 993.147 平方米。利用悬链线模型求解锚链形状，与原结果进行比较，发现图形相近，验证模型正确。对模型进行灵敏度分析，验证模型的鲁棒性较好。

1.1 第一次遍历寻找浮标吃水深度MATLAB程序

%% 浮标1
g=9.807;%重力加速度
p=1025;%海水密度
m1=1000;%浮标质量
v=12;%海面风速count0=0;
HH=zeros(10,1);
for h1=0.50:0.01:0.75 %吃水深度（遍历寻找最优解）
count0=count0+1;%循环计数V1=pi*h1;%浮标吃水体积
syms fT1;%拉力T（下同）
syms fsi1;%角度SITA（下同）G1=m1*g;
B1=p*g*V1;
Ffeng=0.625*((2-h1)*2)*v^2;
eq11=B1-G1-fT1*cosd(fsi1);
eq12=Ffeng-fT1*sind(fsi1);
[fT1,fsi1]=solve(eq11,eq12,fT1,fsi1);
index1=find(fT1>0);
T1=double(fT1(index1));
si1=double(fsi1(index1));
r1=0;%横坐标长度%% 钢管2
m2=10;%每节钢管质量
V2=pi*0.025^2*1;%钢管体积
L2=1;
G2=m2*g;
B2=p*g*V2;T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=(T2(n)*sind(si2(n)))/(sind(si2(n+1)));a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end
% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度
r2=abs(L2*sind(a2));%横坐标长度%% 重物球4
m4=1200;%重物球质量
p4=7850;%钢的密度
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
m3=100;%钢桶质量
V3=pi*0.15^2*1;%钢桶体积
L3=1;
G3=m3*g;
B3=p*g*V3;T31=T2(5);%取第2段的值
si31=si2(5);
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31))^2);
si32=atand((T31*sind(si31))/(T31*cosd(si31)+B3-G3-T41));
syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index3=find(a3>0);
a3=double(a3(index3));h3=L3*cosd(a3);%高度
r3=abs(L3*sind(a3));%横坐标长度%% 各锚链节5
L5=0.105;%单个锚链长度
m5=7*L5;%单个锚链质量
p5=7850;%钢的密度
V5=m5/p5;%单个锚链体积
G5=m5*g;
B5=p*g*V5;T5=zeros(211,1);
si5=zeros(211,1);
a5=zeros(211,1);
h5=zeros(211,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:210si5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end
%% 绘图6
x=zeros(211,1);
y=zeros(211,1);
x(1)=0;
y(1)=0;
count=211;
count00=0;%链接漂浮个数
for ii=1:210if a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(211,1);%总高度
r0=r1+sum(r2)+r3+x(211,1);%游动半径
HH(count0)=H;
end
% plot(x,y,'linewidth',2)
% axis([0 20 0 20])
% xlabel('x','FontSize',28);
% ylabel('y','FontSize',28);
% set(gca,'FontSize',28,'linewidth',1);

1.2 第二次遍历寻找浮标吃水深度MATLAB程序

%% 浮标1
g=9.807;%重力加速度
p=1025;%海水密度
m1=1000;%浮标质量
v=12;%海面风速count0=0;
HH=zeros(10,1);
for h1=0.680:0.001:0.690 %吃水深度（遍历寻找最优解）
count0=count0+1;%循环计数V1=pi*h1;%浮标吃水体积
syms fT1;%拉力T（下同）
syms fsi1;%角度SITA（下同）G1=m1*g;
B1=p*g*V1;
Ffeng=0.625*((2-h1)*2)*v^2;
eq11=B1-G1-fT1*cosd(fsi1);
eq12=Ffeng-fT1*sind(fsi1);
[fT1,fsi1]=solve(eq11,eq12,fT1,fsi1);
index1=find(fT1>0);
T1=double(fT1(index1));
si1=double(fsi1(index1));
r1=0;%横坐标长度%% 钢管2
m2=10;%每节钢管质量
V2=pi*0.025^2*1;%钢管体积
L2=1;
G2=m2*g;
B2=p*g*V2;T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=(T2(n)*sind(si2(n)))/(sind(si2(n+1)));a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end
% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度
r2=abs(L2*sind(a2));%横坐标长度%% 重物球4
m4=1200;%重物球质量
p4=7850;%钢的密度
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
m3=100;%钢桶质量
V3=pi*0.15^2*1;%钢桶体积
L3=1;
G3=m3*g;
B3=p*g*V3;T31=T2(5);%取第2段的值
si31=si2(5);
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31))^2);
si32=atand((T31*sind(si31))/(T31*cosd(si31)+B3-G3-T41));
syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index3=find(a3>0);
a3=double(a3(index3));h3=L3*cosd(a3);%高度
r3=abs(L3*sind(a3));%横坐标长度%% 各锚链节5
L5=0.105;%单个锚链长度
m5=7*L5;%单个锚链质量
p5=7850;%钢的密度
V5=m5/p5;%单个锚链体积
G5=m5*g;
B5=p*g*V5;T5=zeros(211,1);
si5=zeros(211,1);
a5=zeros(211,1);
h5=zeros(211,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:210si5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end
%% 绘图6
x=zeros(211,1);
y=zeros(211,1);
x(1)=0;
y(1)=0;
count=211;
count00=0;%链接漂浮个数
for ii=1:210if a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(211,1);%总高度
r0=r1+sum(r2)+r3+x(211,1);%游动半径
HH(count0)=H;
end
% plot(x,y,'linewidth',2)
% axis([0 20 0 20])
% xlabel('x','FontSize',28);
% ylabel('y','FontSize',28);
% set(gca,'FontSize',28,'linewidth',1);

1.3 代入遍历结果浮标吃水深度MATLAB程序

%% 浮标1
g=9.807;%重力加速度
p=1025;%海水密度
m1=1000;%浮标质量
v=12;%海面风速h1=0.683;%吃水深度V1=pi*h1;%浮标吃水体积
syms fT1;%拉力T（下同）
syms fsi1;%角度SITA（下同）G1=m1*g;
B1=p*g*V1;
Ffeng=0.625*((2-h1)*2)*v^2;
eq11=B1-G1-fT1*cosd(fsi1);
eq12=Ffeng-fT1*sind(fsi1);
[fT1,fsi1]=solve(eq11,eq12,fT1,fsi1);
index1=find(fT1>0);
T1=double(fT1(index1));
si1=double(fsi1(index1));
r1=0;%横坐标长度%% 钢管2
m2=10;%每节钢管质量
V2=pi*0.025^2*1;%钢管体积
L2=1;
G2=m2*g;
B2=p*g*V2;T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=(T2(n)*sind(si2(n)))/(sind(si2(n+1)));a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end
% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度
r2=abs(L2*sind(a2));%横坐标长度%% 重物球4
m4=1200;%重物球质量
p4=7850;%钢的密度
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
m3=100;%钢桶质量
V3=pi*0.15^2*1;%钢桶体积
L3=1;
G3=m3*g;
B3=p*g*V3;T31=T2(5);%取第2段的值
si31=si2(5);
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31))^2);
si32=atand((T31*sind(si31))/(T31*cosd(si31)+B3-G3-T41));
syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index3=find(a3>0);
a3=double(a3(index3));h3=L3*cosd(a3);%高度
r3=abs(L3*sind(a3));%横坐标长度%% 各锚链节5
L5=0.105;%单个锚链长度
m5=7*L5;%单个锚链质量
p5=7850;%钢的密度
V5=m5/p5;%单个锚链体积
G5=m5*g;
B5=p*g*V5;T5=zeros(211,1);
si5=zeros(211,1);
a5=zeros(211,1);
h5=zeros(211,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:210si5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end
%% 绘图6
x=zeros(211,1);
y=zeros(211,1);
x(1)=0;
y(1)=0;
count=211;
count00=0;%链接漂浮个数
for ii=1:210if a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(211,1);%总高度
r0=r1+sum(r2)+r3+x(211,1);%游动半径
plot(x,y,'linewidth',2)
axis([0 20 0 20])
xlabel('x','FontSize',28);
ylabel('y','FontSize',28);
set(gca,'FontSize',28,'linewidth',1);%% 导出数据7
filterName=["钢管角度","钢桶角度","浮标吃水深度","浮标游动半径","锚链末端夹角","链结漂浮个数","系统总高度"];
xlswrite('problem1.xlsx',filterName,'v=12','A1');
xlswrite('problem1.xlsx',a2,'v=12','A2');
xlswrite('problem1.xlsx',a3,'v=12','B2');
xlswrite('problem1.xlsx',h1,'v=12','C2');
xlswrite('problem1.xlsx',r0,'v=12','D2');
xlswrite('problem1.xlsx',90-a5(210),'v=12','E2');
xlswrite('problem1.xlsx',count00,'v=12','F2');
xlswrite('problem1.xlsx',H,'v=12','G2');

二、第2问

针对问题二，首先判断在问题一的假设下，系泊系统对风速为 36m/s 的适 应能力，得出系统不满足题目要求。建立整体力学模型和基于重物球质量的系 泊系统设计优化模型，将钢桶的倾斜角度、浮标的吃水深度、浮标的游动区域、锚链在锚点与海床的夹角、重物球的质量均达到最小的多目标优化转化为 单目标进行优化。根据上下边界条件，求得重物球质量的调节范围，对范围变 步长划分进行循环遍历，寻找其最优值，使系统的稳定性达到最大。结果得到 重物球质量的最优值为 5527kg，系统的稳定性最好。

2.1 遍历寻找满足条件重物球质量最小解MATLAB程序

%% 求极限值
g=9.807;%重力加速度
p=1025;%海水密度
v=36;%海面风速h1=2;%先取最大值就极值
m1=1000;%浮标质量
V=pi*h1;%浮标体积
B=p*g*V;%浮标最大浮力
G=m1*g;%浮标重力
FFfeng=0.625*((2-h1)*2)*v^2;%F风m2=10;%每节钢管质量
V2=pi*0.025^2*1;%钢管体积
L2=1;
G2=m2*g;
B2=p*g*V2;m3=100;%钢桶质量
V3=pi*0.15^2*1;%钢桶体积
L3=1;
G3=m3*g;
B3=p*g*V3;L5=0.105;%单个锚链长度
m5=7*L5;%单个锚链质量
p5=7850;%钢的密度
V5=m5/p5;%单个锚链体积
G5=m5*g;
B5=p*g*V5;syms fm4
p4=7850;%钢的密度
V4=fm4/p4;%重物球体积
G4=fm4*g;
B4=p*g*V4;
T41=G4-B4;
eq0=B+B2*4+B3+B4+B5*210-FFfeng*tand(16)-(G+G2*4+G3+G5*210)-fm4*g;%G铁max=F浮max-F风*tan16°-G其他
[fm4]=solve(eq0,fm4);
fm4=double(fm4);
fm4=round(fm4);
%% 总
for iii=1:200
m4=2200+iii*25; %循环遍历
%% 浮标1
count0=0;
HH=zeros(10,1);
%% 重复遍历1
for h1=0.70:0.1:1.50 %吃水深度（遍历寻找最优解）
count0=count0+1;%循环计数m1=1000;%浮标质量
V1=pi*h1;%浮标体积
B1=p*g*V1;%浮标最大浮力
G1=m1*g;%浮标重力
Ffeng=0.625*((2-h1)*2)*v^2;%F风syms fT1;%拉力T（下同）
syms fsi1;%角度SITA（下同）eq11=B1-G1-fT1*cosd(fsi1);
eq12=Ffeng-fT1*sind(fsi1);
[fT1,fsi1]=solve(eq11,eq12,fT1,fsi1);
index1=find(fT1>0);
T1=double(fT1(index1));
si1=double(fsi1(index1));%% 钢管2
T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=(T2(n)*sind(si2(n)))/(sind(si2(n+1)));a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end
% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度%% 重物球4
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
T31=T2(5);%取第2段的值
si31=si2(5);
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31))^2);
si32=atand((T31*sind(si31))/(T31*cosd(si31)+B3-G3-T41));
syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index31=find(a3>0);
a3=double(a3(index31));h3=L3*cosd(a3);%高度%% 各锚链节5
T5=zeros(211,1);
si5=zeros(211,1);
a5=zeros(211,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:210si5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end%% 绘图6
x=zeros(211,1);
y=zeros(211,1);
x(1)=0;
y(1)=0;
count=211;
count00=0;%链接漂浮个数
for ii=1:210if a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(211,1);%总高度
HH(count0)=H;
% plot(x,y,'linewidth',2)
% axis([0 20 0 20])
% xlabel('x','FontSize',28);
% ylabel('y','FontSize',28);
% set(gca,'FontSize',28,'linewidth',1);
endcount1=10;
for n=1:9count1=count1-1;if HH(count1)<18new1=0.7+0.1*(count1-1);breakend
endcount0=0;
HH=zeros(10,1);
%% 重复遍历2
for h1=new1:0.01:new1+0.1 %吃水深度（遍历寻找最优解）
count0=count0+1;%循环计数m1=1000;%浮标质量
V1=pi*h1;%浮标体积
B1=p*g*V1;%浮标最大浮力
G1=m1*g;%浮标重力
Ffeng=0.625*((2-h1)*2)*v^2;%F风syms fT1;%拉力T（下同）
syms fsi1;%角度SITA（下同）eq11=B1-G1-fT1*cosd(fsi1);
eq12=Ffeng-fT1*sind(fsi1);
[fT1,fsi1]=solve(eq11,eq12,fT1,fsi1);
index1=find(fT1>0);
T1=double(fT1(index1));
si1=double(fsi1(index1));%% 钢管2
T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=(T2(n)*sind(si2(n)))/(sind(si2(n+1)));a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end
% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度%% 重物球4
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
T31=T2(5);%取第2段的值
si31=si2(5);
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31))^2);
si32=atand((T31*sind(si31))/(T31*cosd(si31)+B3-G3-T41));
syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index32=find(a3>0);
a3=double(a3(index32));h3=L3*cosd(a3);%高度%% 各锚链节5
T5=zeros(211,1);
si5=zeros(211,1);
a5=zeros(211,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:210si5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end%% 绘图6
x=zeros(211,1);
y=zeros(211,1);
x(1)=0;
y(1)=0;
count=211;
count00=0;%链接漂浮个数
for ii=1:210if a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(211,1);%总高度
HH(count0)=H;
% plot(x,y,'linewidth',2)
% axis([0 20 0 20])
% xlabel('x','FontSize',28);
% ylabel('y','FontSize',28);
% set(gca,'FontSize',28,'linewidth',1);
endcount2=12;
for n=1:11count2=count2-1;if HH(count2)<18new2=new1+0.01*(count2-1);breakend
endcount0=0;
HH=zeros(10,1);
%% 重复遍历3
for h1=new2:0.001:new2+0.01 %吃水深度（遍历寻找最优解）
count0=count0+1;%循环计数m1=1000;%浮标质量
V1=pi*h1;%浮标体积
B1=p*g*V1;%浮标最大浮力
G1=m1*g;%浮标重力
Ffeng=0.625*((2-h1)*2)*v^2;%F风syms fT1;%拉力T（下同）
syms fsi1;%角度SITA（下同）eq11=B1-G1-fT1*cosd(fsi1);
eq12=Ffeng-fT1*sind(fsi1);
[fT1,fsi1]=solve(eq11,eq12,fT1,fsi1);
index1=find(fT1>0);
T1=double(fT1(index1));
si1=double(fsi1(index1));%% 钢管2
T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=(T2(n)*sind(si2(n)))/(sind(si2(n+1)));a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end
% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度%% 重物球4
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
T31=T2(5);%取第2段的值
si31=si2(5);
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31))^2);
si32=atand((T31*sind(si31))/(T31*cosd(si31)+B3-G3-T41));
syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index33=find(a3>0);
a3=double(a3(index33));h3=L3*cosd(a3);%高度%% 各锚链节5
T5=zeros(211,1);
si5=zeros(211,1);
a5=zeros(211,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:210si5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end
%% 绘图6
x=zeros(211,1);
y=zeros(211,1);
x(1)=0;
y(1)=0;
count=211;
count00=0;%链接漂浮个数
for ii=1:210if a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(211,1);%总高度HH(count0)=H;% plot(x,y,'linewidth',2)
% axis([0 20 0 20])
% xlabel('x','FontSize',28);
% ylabel('y','FontSize',28);
% set(gca,'FontSize',28,'linewidth',1);
end
count3=12;
for n=1:11count3=count3-1;if HH(count3)<18new31=abs(HH(count3)-18);new32=abs(HH(count3+1)-18);if  new31<new32new3= new2+0.001*(count3-1);elsenew3= new2+0.001*(count3);endbreakend
endcount0=0;
HH=zeros(10,1);
%% 重复计算4
h1=new3;%吃水深度
m1=1000;%浮标质量
V1=pi*h1;%浮标体积
B1=p*g*V1;%浮标最大浮力
G1=m1*g;%浮标重力
Ffeng=0.625*((2-h1)*2)*v^2;%F风syms fT1;%拉力T（下同）
syms fsi1;%角度SITA（下同）eq11=B1-G1-fT1*cosd(fsi1);
eq12=Ffeng-fT1*sind(fsi1);
[fT1,fsi1]=solve(eq11,eq12,fT1,fsi1);
index1=find(fT1>0);
T1=double(fT1(index1));
si1=double(fsi1(index1));
r1=0;%横坐标长度%% 钢管2
T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=(T2(n)*sind(si2(n)))/(sind(si2(n+1)));a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end
% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度
r2=abs(L2*sind(a2));%横坐标长度%% 重物球4
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
T31=T2(5);%取第2段的值
si31=si2(5);
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31))^2);
si32=atand((T31*sind(si31))/(T31*cosd(si31)+B3-G3-T41));
syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index3=find(a3>0);
a3=double(a3(index3));h3=L3*cosd(a3);%高度
r3=abs(L3*sind(a3));%横坐标长度%% 各锚链节5
T5=zeros(211,1);
si5=zeros(211,1);
a5=zeros(211,1);
h5=zeros(211,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:210si5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end
%% 绘图6
x=zeros(211,1);
y=zeros(211,1);
x(1)=0;
y(1)=0;
count=211;
count00=0;%链接漂浮个数
for ii=1:210if a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(211,1);%总高度
r0=r1+sum(r2)+r3+x(211,1);%游动半径
%%
if a3<5&&a5(210)>74break
end
end
plot(x,y,'linewidth',2)
axis([0 20 0 20])
xlabel('x','FontSize',28);
ylabel('y','FontSize',28);
set(gca,'FontSize',28,'linewidth',1);%% 导出数据7
filterName=["钢管角度","钢桶角度","浮标吃水深度","浮标游动半径","锚链末端夹角","链结漂浮个数","系统总高度","重物质量"];
xlswrite('problem2（优化前）.xlsx',filterName,'v=36','A1');
xlswrite('problem2（优化前）.xlsx',a2,'v=36','A2');
xlswrite('problem2（优化前）.xlsx',a3,'v=36','B2');
xlswrite('problem2（优化前）.xlsx',h1,'v=36','C2');
xlswrite('problem2（优化前）.xlsx',r0,'v=36','D2');
xlswrite('problem2（优化前）.xlsx',90-a5(210),'v=36','E2');
xlswrite('problem2（优化前）.xlsx',count00,'v=36','F2');
xlswrite('problem2（优化前）.xlsx',H,'v=36','G2');
xlswrite('problem2（优化前）.xlsx',m4,'v=36','H2');

2.2 加入多目标优化遍历寻找满足条件重物球质量MATLAB程序

%% 求极限值
g=9.807;%重力加速度
p=1025;%海水密度
v=36;%海面风速h1=2;%先取最大值就极值
m1=1000;%浮标质量
V=pi*h1;%浮标体积
B=p*g*V;%浮标最大浮力
G=m1*g;%浮标重力
FFfeng=0.625*((2-h1)*2)*v^2;%F风m2=10;%每节钢管质量
V2=pi*0.025^2*1;%钢管体积
L2=1;
G2=m2*g;
B2=p*g*V2;m3=100;%钢桶质量
V3=pi*0.15^2*1;%钢桶体积
L3=1;
G3=m3*g;
B3=p*g*V3;L5=0.105;%单个锚链长度
m5=7*L5;%单个锚链质量
p5=7850;%钢的密度
V5=m5/p5;%单个锚链体积
G5=m5*g;
B5=p*g*V5;syms fm4
p4=7850;%钢的密度
V4=fm4/p4;%重物球体积
G4=fm4*g;
B4=p*g*V4;
T41=G4-B4;
eq0=B+B2*4+B3+B4+B5*210-FFfeng*tand(16)-(G+G2*4+G3+G5*210)-fm4*g;%G铁max=F浮max-F风*tan16°-G其他
[fm4]=solve(eq0,fm4);
fm4=double(fm4);
fm4=round(fm4);
%% 总
ei=0;
mi=zeros(15,1);
for iii=1:15
m4=2225+iii*(fm4-2225)/15; %循环遍历
%% 浮标1
count0=0;
HH=zeros(10,1);
%% 重复遍历1
for h1=0.70:0.1:1.90 %吃水深度（遍历寻找最优解）
count0=count0+1;%循环计数m1=1000;%浮标质量
V1=pi*h1;%浮标体积
B1=p*g*V1;%浮标最大浮力
G1=m1*g;%浮标重力
Ffeng=0.625*((2-h1)*2)*v^2;%F风syms fT1;%拉力T（下同）
syms fsi1;%角度SITA（下同）eq11=B1-G1-fT1*cosd(fsi1);
eq12=Ffeng-fT1*sind(fsi1);
[fT1,fsi1]=solve(eq11,eq12,fT1,fsi1);
index1=find(fT1>0);
T1=double(fT1(index1));
si1=double(fsi1(index1));%% 钢管2
T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=(T2(n)*sind(si2(n)))/(sind(si2(n+1)));a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end
% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度%% 重物球4
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
T31=T2(5);%取第2段的值
si31=si2(5);
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31))^2);
si32=atand((T31*sind(si31))/(T31*cosd(si31)+B3-G3-T41));
syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index31=find(a3>0);
a3=double(a3(index31));h3=L3*cosd(a3);%高度%% 各锚链节5
T5=zeros(211,1);
si5=zeros(211,1);
a5=zeros(211,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:210si5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end%% 绘图6
x=zeros(211,1);
y=zeros(211,1);
x(1)=0;
y(1)=0;
count=211;
count00=0;%链接漂浮个数
for ii=1:210if a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(211,1);%总高度
HH(count0)=H;
% plot(x,y,'linewidth',2)
% axis([0 20 0 20])
% xlabel('x','FontSize',28);
% ylabel('y','FontSize',28);
% set(gca,'FontSize',28,'linewidth',1);
endcount1=14;
for n=1:13count1=count1-1;if HH(count1)<18new1=0.7+0.1*(count1-1);breakend
endcount0=0;
HH=zeros(10,1);
%% 重复遍历2
for h1=new1:0.01:new1+0.1 %吃水深度（遍历寻找最优解）
count0=count0+1;%循环计数m1=1000;%浮标质量
V1=pi*h1;%浮标体积
B1=p*g*V1;%浮标最大浮力
G1=m1*g;%浮标重力
Ffeng=0.625*((2-h1)*2)*v^2;%F风syms fT1;%拉力T（下同）
syms fsi1;%角度SITA（下同）eq11=B1-G1-fT1*cosd(fsi1);
eq12=Ffeng-fT1*sind(fsi1);
[fT1,fsi1]=solve(eq11,eq12,fT1,fsi1);
index1=find(fT1>0);
T1=double(fT1(index1));
si1=double(fsi1(index1));%% 钢管2
T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=(T2(n)*sind(si2(n)))/(sind(si2(n+1)));a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end
% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度%% 重物球4
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
T31=T2(5);%取第2段的值
si31=si2(5);
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31))^2);
si32=atand((T31*sind(si31))/(T31*cosd(si31)+B3-G3-T41));
syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index32=find(a3>0);
a3=double(a3(index32));h3=L3*cosd(a3);%高度%% 各锚链节5
T5=zeros(211,1);
si5=zeros(211,1);
a5=zeros(211,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:210si5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end%% 绘图6
x=zeros(211,1);
y=zeros(211,1);
x(1)=0;
y(1)=0;
count=211;
count00=0;%链接漂浮个数
for ii=1:210if a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(211,1);%总高度
HH(count0)=H;
% plot(x,y,'linewidth',2)
% axis([0 20 0 20])
% xlabel('x','FontSize',28);
% ylabel('y','FontSize',28);
% set(gca,'FontSize',28,'linewidth',1);
endcount2=12;
for n=1:11count2=count2-1;if HH(count2)<18new2=new1+0.01*(count2-1);breakend
endcount0=0;
HH=zeros(10,1);
%% 重复遍历3
for h1=new2:0.001:new2+0.01 %吃水深度（遍历寻找最优解）
count0=count0+1;%循环计数m1=1000;%浮标质量
V1=pi*h1;%浮标体积
B1=p*g*V1;%浮标最大浮力
G1=m1*g;%浮标重力
Ffeng=0.625*((2-h1)*2)*v^2;%F风syms fT1;%拉力T（下同）
syms fsi1;%角度SITA（下同）eq11=B1-G1-fT1*cosd(fsi1);
eq12=Ffeng-fT1*sind(fsi1);
[fT1,fsi1]=solve(eq11,eq12,fT1,fsi1);
index1=find(fT1>0);
T1=double(fT1(index1));
si1=double(fsi1(index1));%% 钢管2
T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=(T2(n)*sind(si2(n)))/(sind(si2(n+1)));a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end
% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度%% 重物球4
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
T31=T2(5);%取第2段的值
si31=si2(5);
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31))^2);
si32=atand((T31*sind(si31))/(T31*cosd(si31)+B3-G3-T41));
syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index33=find(a3>0);
a3=double(a3(index33));h3=L3*cosd(a3);%高度%% 各锚链节5
T5=zeros(211,1);
si5=zeros(211,1);
a5=zeros(211,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:210si5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end
%% 绘图6
x=zeros(211,1);
y=zeros(211,1);
x(1)=0;
y(1)=0;
count=211;
count00=0;%链接漂浮个数
for ii=1:210if a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(211,1);%总高度HH(count0)=H;% plot(x,y,'linewidth',2)
% axis([0 20 0 20])
% xlabel('x','FontSize',28);
% ylabel('y','FontSize',28);
% set(gca,'FontSize',28,'linewidth',1);
end
count3=12;
for n=1:11count3=count3-1;if HH(count3)<18if count3==11new3= new2+0.001*(count3-1);elsenew31=abs(HH(count3)-18);new32=abs(HH(count3+1)-18);if   new31<new32new3= new2+0.001*(count3-1);elsenew3= new2+0.001*(count3);endendbreakend
endcount0=0;
HH=zeros(10,1);
%% 重复计算4
h1=new3;%吃水深度
m1=1000;%浮标质量
V1=pi*h1;%浮标体积
B1=p*g*V1;%浮标最大浮力
G1=m1*g;%浮标重力
Ffeng=0.625*((2-h1)*2)*v^2;%F风syms fT1;%拉力T（下同）
syms fsi1;%角度SITA（下同）eq11=B1-G1-fT1*cosd(fsi1);
eq12=Ffeng-fT1*sind(fsi1);
[fT1,fsi1]=solve(eq11,eq12,fT1,fsi1);
index1=find(fT1>0);
T1=double(fT1(index1));
si1=double(fsi1(index1));
r1=0;%横坐标长度%% 钢管2
T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=(T2(n)*sind(si2(n)))/(sind(si2(n+1)));a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end
% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度
r2=abs(L2*sind(a2));%横坐标长度%% 重物球4
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
T31=T2(5);%取第2段的值
si31=si2(5);
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31))^2);
si32=atand((T31*sind(si31))/(T31*cosd(si31)+B3-G3-T41));
syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index3=find(a3>0);
a3=double(a3(index3));h3=L3*cosd(a3);%高度
r3=abs(L3*sind(a3));%横坐标长度%% 各锚链节5
T5=zeros(211,1);
si5=zeros(211,1);
a5=zeros(211,1);
h5=zeros(211,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:210si5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end
%% 绘图6
x=zeros(211,1);
y=zeros(211,1);
x(1)=0;
y(1)=0;
count=211;
count00=0;%链接漂浮个数
for ii=1:210if a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(211,1);%总高度
r0=r1+sum(r2)+r3+x(211,1);%游动半径
%% 多目标优化
ei=ei+1;
a3max=5;%角度a3极限值
h1max=2;%吃水深度h1极限值
r0max=27.05;%R0游动半径极限值
a5210max=90-16;%锚链与锚夹角a5(210)极限值
mi(ei)=a3/a3max+h1/h1max+r0/r0max+(90-a5(210))/(a5210max)+m4/(m4+fm4);
end
new4=find(mi==min(mi));
new4=2225+new4*(fm4-2225)/15;

2.3 代入遍历优化结果重物球质量MATLAB程序

%% 求极限值
g=9.807;%重力加速度
p=1025;%海水密度
v=36;%海面风速h1=2;%先取最大值就极值
m1=1000;%浮标质量
V=pi*h1;%浮标体积
B=p*g*V;%浮标最大浮力
G=m1*g;%浮标重力
FFfeng=0.625*((2-h1)*2)*v^2;%F风m2=10;%每节钢管质量
V2=pi*0.025^2*1;%钢管体积
L2=1;
G2=m2*g;
B2=p*g*V2;m3=100;%钢桶质量
V3=pi*0.15^2*1;%钢桶体积
L3=1;
G3=m3*g;
B3=p*g*V3;L5=0.105;%单个锚链长度
m5=7*L5;%单个锚链质量
p5=7850;%钢的密度
V5=m5/p5;%单个锚链体积
G5=m5*g;
B5=p*g*V5;syms fm4
p4=7850;%钢的密度
V4=fm4/p4;%重物球体积
G4=fm4*g;
B4=p*g*V4;
T41=G4-B4;
eq0=B+B2*4+B3+B4+B5*210-FFfeng*tand(16)-(G+G2*4+G3+G5*210)-fm4*g;%G铁max=F浮max-F风*tan16°-G其他
[fm4]=solve(eq0,fm4);
fm4=double(fm4);
fm4=round(fm4);
%% 总
m4=5527; %最优解
%% 浮标1
count0=0;
HH=zeros(10,1);
%% 重复遍历1
for h1=0.70:0.1:1.90 %吃水深度
count0=count0+1;%循环计数m1=1000;%浮标质量
V1=pi*h1;%浮标体积
B1=p*g*V1;%浮标最大浮力
G1=m1*g;%浮标重力
Ffeng=0.625*((2-h1)*2)*v^2;%F风syms fT1;%拉力T（下同）
syms fsi1;%角度SITA（下同）eq11=B1-G1-fT1*cosd(fsi1);
eq12=Ffeng-fT1*sind(fsi1);
[fT1,fsi1]=solve(eq11,eq12,fT1,fsi1);
index1=find(fT1>0);
T1=double(fT1(index1));
si1=double(fsi1(index1));%% 钢管2
T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=(T2(n)*sind(si2(n)))/(sind(si2(n+1)));a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end
% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度%% 重物球4
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
T31=T2(5);%取第2段的值
si31=si2(5);
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31))^2);
si32=atand((T31*sind(si31))/(T31*cosd(si31)+B3-G3-T41));
syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index31=find(a3>0);
a3=double(a3(index31));h3=L3*cosd(a3);%高度%% 各锚链节5
T5=zeros(211,1);
si5=zeros(211,1);
a5=zeros(211,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:210si5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end%% 绘图6
x=zeros(211,1);
y=zeros(211,1);
x(1)=0;
y(1)=0;
count=211;
count00=0;%链接漂浮个数
for ii=1:210if a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(211,1);%总高度
HH(count0)=H;
% plot(x,y,'linewidth',2)
% axis([0 20 0 20])
% xlabel('x','FontSize',28);
% ylabel('y','FontSize',28);
% set(gca,'FontSize',28,'linewidth',1);
endcount1=14;
for n=1:13count1=count1-1;if HH(count1)<18new1=0.7+0.1*(count1-1);breakend
endcount0=0;
HH=zeros(10,1);
%% 重复遍历2
for h1=new1:0.01:new1+0.1 %吃水深度（遍历寻找最优解）
count0=count0+1;%循环计数m1=1000;%浮标质量
V1=pi*h1;%浮标体积
B1=p*g*V1;%浮标最大浮力
G1=m1*g;%浮标重力
Ffeng=0.625*((2-h1)*2)*v^2;%F风syms fT1;%拉力T（下同）
syms fsi1;%角度SITA（下同）eq11=B1-G1-fT1*cosd(fsi1);
eq12=Ffeng-fT1*sind(fsi1);
[fT1,fsi1]=solve(eq11,eq12,fT1,fsi1);
index1=find(fT1>0);
T1=double(fT1(index1));
si1=double(fsi1(index1));%% 钢管2
T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=(T2(n)*sind(si2(n)))/(sind(si2(n+1)));a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end
% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度%% 重物球4
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
T31=T2(5);%取第2段的值
si31=si2(5);
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31))^2);
si32=atand((T31*sind(si31))/(T31*cosd(si31)+B3-G3-T41));
syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index32=find(a3>0);
a3=double(a3(index32));h3=L3*cosd(a3);%高度%% 各锚链节5
T5=zeros(211,1);
si5=zeros(211,1);
a5=zeros(211,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:210si5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end%% 绘图6
x=zeros(211,1);
y=zeros(211,1);
x(1)=0;
y(1)=0;
count=211;
count00=0;%链接漂浮个数
for ii=1:210if a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(211,1);%总高度
HH(count0)=H;
% plot(x,y,'linewidth',2)
% axis([0 20 0 20])
% xlabel('x','FontSize',28);
% ylabel('y','FontSize',28);
% set(gca,'FontSize',28,'linewidth',1);
endcount2=12;
for n=1:11count2=count2-1;if HH(count2)<18new2=new1+0.01*(count2-1);breakend
endcount0=0;
HH=zeros(10,1);
%% 重复遍历3
for h1=new2:0.001:new2+0.01 %吃水深度（遍历寻找最优解）
count0=count0+1;%循环计数m1=1000;%浮标质量
V1=pi*h1;%浮标体积
B1=p*g*V1;%浮标最大浮力
G1=m1*g;%浮标重力
Ffeng=0.625*((2-h1)*2)*v^2;%F风syms fT1;%拉力T（下同）
syms fsi1;%角度SITA（下同）eq11=B1-G1-fT1*cosd(fsi1);
eq12=Ffeng-fT1*sind(fsi1);
[fT1,fsi1]=solve(eq11,eq12,fT1,fsi1);
index1=find(fT1>0);
T1=double(fT1(index1));
si1=double(fsi1(index1));%% 钢管2
T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=(T2(n)*sind(si2(n)))/(sind(si2(n+1)));a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end
% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度%% 重物球4
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
T31=T2(5);%取第2段的值
si31=si2(5);
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31))^2);
si32=atand((T31*sind(si31))/(T31*cosd(si31)+B3-G3-T41));
syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index33=find(a3>0);
a3=double(a3(index33));h3=L3*cosd(a3);%高度%% 各锚链节5
T5=zeros(211,1);
si5=zeros(211,1);
a5=zeros(211,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:210si5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end
%% 绘图6
x=zeros(211,1);
y=zeros(211,1);
x(1)=0;
y(1)=0;
count=211;
count00=0;%链接漂浮个数
for ii=1:210if a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(211,1);%总高度HH(count0)=H;% plot(x,y,'linewidth',2)
% axis([0 20 0 20])
% xlabel('x','FontSize',28);
% ylabel('y','FontSize',28);
% set(gca,'FontSize',28,'linewidth',1);
end
count3=12;
for n=1:11count3=count3-1;if HH(count3)<18if count3==11new3= new2+0.001*(count3);elsenew31=abs(HH(count3)-18);new32=abs(HH(count3+1)-18);if   new31<new32new3= new2+0.001*(count3-1);elsenew3= new2+0.001*(count3);endendbreakend
endcount0=0;
HH=zeros(10,1);
%% 重复计算4
h1=new3;%吃水深度
m1=1000;%浮标质量
V1=pi*h1;%浮标体积
B1=p*g*V1;%浮标最大浮力
G1=m1*g;%浮标重力
Ffeng=0.625*((2-h1)*2)*v^2;%F风syms fT1;%拉力T（下同）
syms fsi1;%角度SITA（下同）eq11=B1-G1-fT1*cosd(fsi1);
eq12=Ffeng-fT1*sind(fsi1);
[fT1,fsi1]=solve(eq11,eq12,fT1,fsi1);
index1=find(fT1>0);
T1=double(fT1(index1));
si1=double(fsi1(index1));
r1=0;%横坐标长度%% 钢管2
T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=(T2(n)*sind(si2(n)))/(sind(si2(n+1)));a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end
% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度
r2=abs(L2*sind(a2));%横坐标长度%% 重物球4
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
T31=T2(5);%取第2段的值
si31=si2(5);
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31))^2);
si32=atand((T31*sind(si31))/(T31*cosd(si31)+B3-G3-T41));
syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index3=find(a3>0);
a3=double(a3(index3));h3=L3*cosd(a3);%高度
r3=abs(L3*sind(a3));%横坐标长度%% 各锚链节5
T5=zeros(211,1);
si5=zeros(211,1);
a5=zeros(211,1);
h5=zeros(211,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:210si5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end
%% 绘图6
x=zeros(211,1);
y=zeros(211,1);
x(1)=0;
y(1)=0;
count=211;
count00=0;%链接漂浮个数
for ii=1:210if a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(211,1);%总高度
r0=r1+sum(r2)+r3+x(211,1);%游动半径plot(x,y,'linewidth',2)
axis([0 20 0 20])
xlabel('x','FontSize',28);
ylabel('y','FontSize',28);
set(gca,'FontSize',28,'linewidth',1);
%% 导出数据7
filterName=["钢管角度","钢桶角度","浮标吃水深度","浮标游动半径","锚链末端夹角","链结漂浮个数","系统总高度","重物质量"];
xlswrite('problem2（优化后）.xlsx',filterName,'v=36','A1');
xlswrite('problem2（优化后）.xlsx',a2,'v=36','A2');
xlswrite('problem2（优化后）.xlsx',a3,'v=36','B2');
xlswrite('problem2（优化后）.xlsx',h1,'v=36','C2');
xlswrite('problem2（优化后）.xlsx',r0,'v=36','D2');
xlswrite('problem2（优化后）.xlsx',90-a5(210),'v=36','E2');
xlswrite('problem2（优化后）.xlsx',count00,'v=36','F2');
xlswrite('problem2（优化后）.xlsx',H,'v=36','G2');
xlswrite('problem2（优化后）.xlsx',m4,'v=36','H2');

三、第3问

针对问题三，增加水流力这一影响因素，对系泊系统各部分受力分析进行调整。利用问题一的方法，选取不同锚链型号，长度以及重物球的质量，进行变步长遍历。沿用问题二建立的优化模型，求解最佳设计方案，使系泊系统能 够适应极端环境。结果得到锚链型号为Ⅴ号，链环个数为 125 个，重物球质量为5000kg，系统的稳定性最好。选取两组不同的水深、海水速度、风速，对设计的系泊系统进行检验，得到系泊系统的稳定性好。

3.1 遍历不同锚链型号与长度以及重物球的质量MATLAB程序

%% 初始化
g=9.807;%重力加速度
p=1025;%海水密度
v=36;%海面风速
v0=1.5;%水速h1=2;%先取最大值就极值
m1=1000;%浮标质量
V=pi*h1;%浮标体积
B=p*g*V;%浮标最大浮力
G=m1*g;%浮标重力
FFfeng=0.625*((2-h1)*2)*v^2;%F风m2=10;%每节钢管质量
V2=pi*0.025^2*1;%钢管体积
L2=1;
G2=m2*g;
B2=p*g*V2;m3=100;%钢桶质量
V3=pi*0.15^2*1;%钢桶体积
L3=1;
G3=m3*g;
B3=p*g*V3;L5=0.105;%单个锚链长度
m5=7*L5;%单个锚链质量
p5=7850;%钢的密度
V5=m5/p5;%单个锚链体积
G5=m5*g;
B5=p*g*V5;N=[0.078,3.2;0.105,7;0.120,12.5;0.150,19.5;0.180,28.12];
mi=zeros(10,1);
ei=0;
for nn=1:5
L5=N(nn,1);%单个锚链长度
m5=N(nn,2)*L5;%单个锚链质量
p5=7850;%钢的密度
V5=m5/p5;%单个锚链体积
G5=m5*g;
B5=p*g*V5;
S=(20-5)/L5;
S=round(S);
for s=S:25:150+S
syms ffm4
p4=7850;%钢的密度
V4=ffm4/p4;%重物球体积
G4=ffm4*g;
B4=p*g*V4;
T41=G4-B4;
eq0=B+B2*4+B3+B4+B5*S-FFfeng*tand(16)-(G+G2*4+G3+G5*S)-ffm4*g;%G铁max=F浮max-F风*tan16°-G其他
[ffm4]=solve(eq0,ffm4);
fm4=double(ffm4);
fm4=round(fm4);
for iii=1:15
m4=2225+iii*(fm4-2225)/15; %循环遍历
%% 浮标1
count0=0;
HH=zeros(10,1);
%% 重复遍历1
for h1=0.70:0.1:1.90 %吃水深度（遍历寻找最优解）
count0=count0+1;%循环计数m1=1000;%浮标质量
V1=pi*h1;%浮标体积
B1=p*g*V1;%浮标最大浮力
G1=m1*g;%浮标重力
S1=2*h1;%受水投影面积
Ffeng=0.625*((2-h1)*2)*v^2;%F风
Fshui1=374*S1*v0^2;si1=atand((Fshui1+Ffeng)/(B1-G1));
T1=sqrt((B1-G1)^2+(Fshui1+Ffeng)^2);% syms fT1;%拉力T（下同）
% syms fsi1;%角度SITA（下同）
% eq11=B1-G1-fT1*cosd(fsi1);
% eq12=Ffeng-fT1*sind(fsi1);
% [fT1,fsi1]=solve(eq11,eq12,fT1,fsi1);
% index1=find(fT1>0);
% T1=double(fT1(index1));
% si1=double(fsi1(index1));%% 钢管2
T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
S2=0.05*1;
Fshui2=374*S2*v0^2;
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n))+Fshui2)/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=sqrt((T2(n)*cosd(si2(n))+B2-G2)^2+(T2(n)*sind(si2(n))+Fshui2)^2);a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end
% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度%% 重物球4
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
T31=T2(5);%取第2段的值
si31=si2(5);
S3=0.3*1;
Fshui3=374*S3*v0^2;
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31)+Fshui3)^2);
si32=atand((T31*sind(si31)+Fshui3)/(T31*cosd(si31)+B3-G3-T41));
syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index31=find(a3>0);
a3=double(a3(index31));h3=L3*cosd(a3);%高度%% 各锚链节5
T5=zeros(S+1,1);
si5=zeros(S+1,1);
a5=zeros(S+1,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:Ssi5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end%% 绘图6
x=zeros(S+1,1);
y=zeros(S+1,1);
x(1)=0;
y(1)=0;
count=S+1;
count00=0;%链接漂浮个数
for ii=1:Sif a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(S+1,1);%总高度
HH(count0)=H;
% plot(x,y,'linewidth',2)
% axis([0 20 0 20])
% xlabel('x','FontSize',28);
% ylabel('y','FontSize',28);
% set(gca,'FontSize',28,'linewidth',1);
end
if min(HH)>20continue
endcount1=14;
for n=1:13count1=count1-1;if HH(count1)<20new1=0.7+0.1*(count1-1);breakend
endcount0=0;
HH=zeros(10,1);
%% 重复遍历2
for h1=new1:0.01:new1+0.1 %吃水深度（遍历寻找最优解）
count0=count0+1;%循环计数m1=1000;%浮标质量
V1=pi*h1;%浮标体积
B1=p*g*V1;%浮标最大浮力
G1=m1*g;%浮标重力
S1=2*h1;%受水投影面积
Ffeng=0.625*((2-h1)*2)*v^2;%F风
Fshui1=374*S1*v0^2;si1=atand((Fshui1+Ffeng)/(B1-G1));
T1=sqrt((B1-G1)^2+(Fshui1+Ffeng)^2);% syms fT1;%拉力T（下同）
% syms fsi1;%角度SITA（下同）
% eq11=B1-G1-fT1*cosd(fsi1);
% eq12=Ffeng-fT1*sind(fsi1);
% [fT1,fsi1]=solve(eq11,eq12,fT1,fsi1);
% index1=find(fT1>0);
% T1=double(fT1(index1));
% si1=double(fsi1(index1));%% 钢管2
T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
S2=0.05*1;
Fshui2=374*S2*v0^2;
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n))+Fshui2)/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=sqrt((T2(n)*cosd(si2(n))+B2-G2)^2+(T2(n)*sind(si2(n))+Fshui2)^2);a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end
% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度%% 重物球4
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
T31=T2(5);%取第2段的值
si31=si2(5);
S3=0.3*1;
Fshui3=374*S3*v0^2;
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31)+Fshui3)^2);
si32=atand((T31*sind(si31)+Fshui3)/(T31*cosd(si31)+B3-G3-T41));
syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index31=find(a3>0);
a3=double(a3(index31));h3=L3*cosd(a3);%高度%% 各锚链节5
T5=zeros(S+1,1);
si5=zeros(S+1,1);
a5=zeros(S+1,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:Ssi5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end%% 绘图6
x=zeros(S+1,1);
y=zeros(S+1,1);
x(1)=0;
y(1)=0;
count=S+1;
count00=0;%链接漂浮个数
for ii=1:Sif a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(S+1,1);%总高度
HH(count0)=H;
% plot(x,y,'linewidth',2)
% axis([0 20 0 20])
% xlabel('x','FontSize',28);
% ylabel('y','FontSize',28);
% set(gca,'FontSize',28,'linewidth',1);
end
if min(HH)>20continue
end
count2=12;
for n=1:11count2=count2-1;if HH(count2)<20new2=new1+0.01*(count2-1);breakend
end
count0=0;
HH=zeros(10,1);
%% 重复遍历3
for h1=new2:0.001:new2+0.01 %吃水深度（遍历寻找最优解）
count0=count0+1;%循环计数m1=1000;%浮标质量
V1=pi*h1;%浮标体积
B1=p*g*V1;%浮标最大浮力
G1=m1*g;%浮标重力
S1=2*h1;%受水投影面积
Ffeng=0.625*((2-h1)*2)*v^2;%F风
Fshui1=374*S1*v0^2;si1=atand((Fshui1+Ffeng)/(B1-G1));
T1=sqrt((B1-G1)^2+(Fshui1+Ffeng)^2);% syms fT1;%拉力T（下同）
% syms fsi1;%角度SITA（下同）
% eq11=B1-G1-fT1*cosd(fsi1);
% eq12=Ffeng-fT1*sind(fsi1);
% [fT1,fsi1]=solve(eq11,eq12,fT1,fsi1);
% index1=find(fT1>0);
% T1=double(fT1(index1));
% si1=double(fsi1(index1));%% 钢管2
T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
S2=0.05*1;
Fshui2=374*S2*v0^2;
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n))+Fshui2)/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=sqrt((T2(n)*cosd(si2(n))+B2-G2)^2+(T2(n)*sind(si2(n))+Fshui2)^2);a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end
% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度%% 重物球4
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
T31=T2(5);%取第2段的值
si31=si2(5);
S3=0.3*1;
Fshui3=374*S3*v0^2;
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31)+Fshui3)^2);
si32=atand((T31*sind(si31)+Fshui3)/(T31*cosd(si31)+B3-G3-T41));
syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index31=a3>0;
a3=double(a3(index31));h3=L3*cosd(a3);%高度%% 各锚链节5
T5=zeros(S+1,1);
si5=zeros(S+1,1);
a5=zeros(S+1,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:Ssi5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end%% 绘图6
x=zeros(S+1,1);
y=zeros(S+1,1);
x(1)=0;
y(1)=0;
count=S+1;
count00=0;%链接漂浮个数
for ii=1:Sif a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(S+1,1);%总高度
HH(count0)=H;
% plot(x,y,'linewidth',2)
% axis([0 20 0 20])
% xlabel('x','FontSize',28);
% ylabel('y','FontSize',28);
% set(gca,'FontSize',28,'linewidth',1);
end
count3=12;
for n=1:11count3=count3-1;if HH(count3)<20if count3==11new3= new2+0.001*(count3-1);elsenew31=abs(HH(count3)-20);new32=abs(HH(count3+1)-20);if   new31<new32new3= new2+0.001*(count3-1);elsenew3= new2+0.001*(count3);endendbreakend
endcount0=0;
HH=zeros(10,1);
%% 重复计算4
h1=new3;
m1=1000;%浮标质量
V1=pi*h1;%浮标体积
B1=p*g*V1;%浮标最大浮力
G1=m1*g;%浮标重力
S1=2*h1;%受水投影面积
Ffeng=0.625*((2-h1)*2)*v^2;%F风
Fshui1=374*S1*v0^2;si1=atand((Fshui1+Ffeng)/(B1-G1));
T1=sqrt((B1-G1)^2+(Fshui1+Ffeng)^2);% syms fT1;%拉力T（下同）
% syms fsi1;%角度SITA（下同）
% eq11=B1-G1-fT1*cosd(fsi1);
% eq12=Ffeng-fT1*sind(fsi1);
% [fT1,fsi1]=solve(eq11,eq12,fT1,fsi1);
% index1=find(fT1>0);
% T1=double(fT1(index1));
% si1=double(fsi1(index1));
r1=0;%横坐标长度
%% 钢管2
T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
S2=0.05*1;
Fshui2=374*S2*v0^2;
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n))+Fshui2)/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=sqrt((T2(n)*cosd(si2(n))+B2-G2)^2+(T2(n)*sind(si2(n))+Fshui2)^2);a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end
% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度
r2=abs(L2*sind(a2));%横坐标长度
%% 重物球4
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
T31=T2(5);%取第2段的值
si31=si2(5);
S3=0.3*1;
Fshui3=374*S3*v0^2;
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31)+Fshui3)^2);
si32=atand((T31*sind(si31)+Fshui3)/(T31*cosd(si31)+B3-G3-T41));
syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index31=find(a3>0);
a3=double(a3(index31));h3=L3*cosd(a3);%高度
r3=abs(L3*sind(a3));%横坐标长度
%% 各锚链节5
T5=zeros(S+1,1);
si5=zeros(S+1,1);
a5=zeros(S+1,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:Ssi5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end%% 绘图计算值6
x=zeros(S+1,1);
y=zeros(S+1,1);
x(1)=0;
y(1)=0;
count=S+1;
count00=0;%链接漂浮个数
for ii=1:Sif a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(S+1,1);%总高度
r0=r1+sum(r2)+r3+x(S+1,1);%游动半径
%% 多目标优化
ei=ei+1;
a3max=5;%角度a3极限值
h1max=2;%吃水深度h1极限值
r0max=27.05;%R0游动半径极限值
a5210max=90-16;%锚链与锚夹角a5(S)极限值
if a3<5&&a5(S)>74mi(ei)=a3/a3max+h1/h1max+r0/r0max+(90-a5(S))/(a5210max)+m4/(m4+fm4);
elsemi(ei)=inf;
end
end
end
end
new4=find(mi==min(mi));

3.2 三次遍历代入锚链型号与长度以及重物球的质量求浮标吃水深度MATLAB程序

3.2.1 第一次遍历

%% 浮标1
g=9.807;%重力加速度
p=1025;%海水密度
m1=1000;%浮标质量
v=36;%海面风速
v0=1.5;%海水水速
nnn=125;%锚链节个数count0=0;
HH=zeros(10,1);
for h1=0.50:0.1:1.9 %吃水深度（遍历寻找最优解）
count0=count0+1;%循环计数m1=1000;%浮标质量
V1=pi*h1;%浮标体积
B1=p*g*V1;%浮标最大浮力
G1=m1*g;%浮标重力
S1=2*h1;%受水投影面积
Ffeng=0.625*((2-h1)*2)*v^2;%F风
Fshui1=374*S1*v0^2;si1=atand((Fshui1+Ffeng)/(B1-G1));
T1=sqrt((B1-G1)^2+(Fshui1+Ffeng)^2);r1=0;%横坐标长度%% 钢管2
m2=10;%每节钢管质量
V2=pi*0.025^2*1;%钢管体积
L2=1;
G2=m2*g;
B2=p*g*V2;T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
S2=0.05*1;
Fshui2=374*S2*v0^2;
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n))+Fshui2)/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=sqrt((T2(n)*cosd(si2(n))+B2-G2)^2+(T2(n)*sind(si2(n))+Fshui2)^2);a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度
r2=abs(L2*sind(a2));%横坐标长度%% 重物球4
m4=5000;%重物球质量
p4=7850;%钢的密度
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
m3=100;%钢桶质量
V3=pi*0.15^2*1;%钢桶体积
L3=1;
G3=m3*g;
B3=p*g*V3;T31=T2(5);%取第2段的值
si31=si2(5);
S3=0.3*1;
Fshui3=374*S3*v0^2;
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31)+Fshui3)^2);
si32=atand((T31*sind(si31)+Fshui3)/(T31*cosd(si31)+B3-G3-T41));syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index3=find(a3>0);
a3=double(a3(index3));h3=L3*cosd(a3);%高度
r3=abs(L3*sind(a3));%横坐标长度%% 各锚链节5
L5=0.180;%单个锚链长度
m5=28.12*L5;%单个锚链质量
p5=7850;%钢的密度
V5=m5/p5;%单个锚链体积
G5=m5*g;
B5=p*g*V5;T5=zeros(nnn+1,1);
si5=zeros(nnn+1,1);
a5=zeros(nnn+1,1);
h5=zeros(nnn+1,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:nnnsi5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end
%% 绘图6
x=zeros(nnn+1,1);
y=zeros(nnn+1,1);
x(1)=0;
y(1)=0;
count=nnn+1;
count00=0;%链接漂浮个数
for ii=1:nnnif a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(nnn+1,1);%总高度
r0=r1+sum(r2)+r3+x(nnn+1,1);%游动半径
HH(count0)=H;
end
% plot(x,y,'linewidth',2)
% axis([0 20 0 20])
% xlabel('x','FontSize',28);
% ylabel('y','FontSize',28);
% set(gca,'FontSize',28,'linewidth',1);

3.2.2 第二次遍历

%% 浮标1
g=9.807;%重力加速度
p=1025;%海水密度
m1=1000;%浮标质量
v=36;%海面风速
v0=1.5;%海水水速
nnn=125;%锚链节个数count0=0;
HH=zeros(10,1);
for h1=1.8:0.01:1.9 %吃水深度（遍历寻找最优解）
count0=count0+1;%循环计数m1=1000;%浮标质量
V1=pi*h1;%浮标体积
B1=p*g*V1;%浮标最大浮力
G1=m1*g;%浮标重力
S1=2*h1;%受水投影面积
Ffeng=0.625*((2-h1)*2)*v^2;%F风
Fshui1=374*S1*v0^2;si1=atand((Fshui1+Ffeng)/(B1-G1));
T1=sqrt((B1-G1)^2+(Fshui1+Ffeng)^2);r1=0;%横坐标长度%% 钢管2
m2=10;%每节钢管质量
V2=pi*0.025^2*1;%钢管体积
L2=1;
G2=m2*g;
B2=p*g*V2;T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
S2=0.05*1;
Fshui2=374*S2*v0^2;
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n))+Fshui2)/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=sqrt((T2(n)*cosd(si2(n))+B2-G2)^2+(T2(n)*sind(si2(n))+Fshui2)^2);a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度
r2=abs(L2*sind(a2));%横坐标长度%% 重物球4
m4=5000;%重物球质量
p4=7850;%钢的密度
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
m3=100;%钢桶质量
V3=pi*0.15^2*1;%钢桶体积
L3=1;
G3=m3*g;
B3=p*g*V3;T31=T2(5);%取第2段的值
si31=si2(5);
S3=0.3*1;
Fshui3=374*S3*v0^2;
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31)+Fshui3)^2);
si32=atand((T31*sind(si31)+Fshui3)/(T31*cosd(si31)+B3-G3-T41));syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index3=find(a3>0);
a3=double(a3(index3));h3=L3*cosd(a3);%高度
r3=abs(L3*sind(a3));%横坐标长度%% 各锚链节5
L5=0.180;%单个锚链长度
m5=28.12*L5;%单个锚链质量
p5=7850;%钢的密度
V5=m5/p5;%单个锚链体积
G5=m5*g;
B5=p*g*V5;T5=zeros(nnn+1,1);
si5=zeros(nnn+1,1);
a5=zeros(nnn+1,1);
h5=zeros(nnn+1,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:nnnsi5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end
%% 绘图6
x=zeros(nnn+1,1);
y=zeros(nnn+1,1);
x(1)=0;
y(1)=0;
count=nnn+1;
count00=0;%链接漂浮个数
for ii=1:nnnif a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(nnn+1,1);%总高度
r0=r1+sum(r2)+r3+x(nnn+1,1);%游动半径
HH(count0)=H;
end
% plot(x,y,'linewidth',2)
% axis([0 20 0 20])
% xlabel('x','FontSize',28);
% ylabel('y','FontSize',28);
% set(gca,'FontSize',28,'linewidth',1);

3.2.3 第三次遍历

%% 浮标1
g=9.807;%重力加速度
p=1025;%海水密度
m1=1000;%浮标质量
v=36;%海面风速
v0=1.5;%海水水速
nnn=125;%锚链节个数count0=0;
HH=zeros(10,1);
for h1=1.86:0.001:1.87 %吃水深度（遍历寻找最优解）
count0=count0+1;%循环计数m1=1000;%浮标质量
V1=pi*h1;%浮标体积
B1=p*g*V1;%浮标最大浮力
G1=m1*g;%浮标重力
S1=2*h1;%受水投影面积
Ffeng=0.625*((2-h1)*2)*v^2;%F风
Fshui1=374*S1*v0^2;si1=atand((Fshui1+Ffeng)/(B1-G1));
T1=sqrt((B1-G1)^2+(Fshui1+Ffeng)^2);r1=0;%横坐标长度%% 钢管2
m2=10;%每节钢管质量
V2=pi*0.025^2*1;%钢管体积
L2=1;
G2=m2*g;
B2=p*g*V2;T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
S2=0.05*1;
Fshui2=374*S2*v0^2;
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n))+Fshui2)/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=sqrt((T2(n)*cosd(si2(n))+B2-G2)^2+(T2(n)*sind(si2(n))+Fshui2)^2);a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度
r2=abs(L2*sind(a2));%横坐标长度%% 重物球4
m4=5000;%重物球质量
p4=7850;%钢的密度
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
m3=100;%钢桶质量
V3=pi*0.15^2*1;%钢桶体积
L3=1;
G3=m3*g;
B3=p*g*V3;T31=T2(5);%取第2段的值
si31=si2(5);
S3=0.3*1;
Fshui3=374*S3*v0^2;
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31)+Fshui3)^2);
si32=atand((T31*sind(si31)+Fshui3)/(T31*cosd(si31)+B3-G3-T41));syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index3=find(a3>0);
a3=double(a3(index3));h3=L3*cosd(a3);%高度
r3=abs(L3*sind(a3));%横坐标长度%% 各锚链节5
L5=0.180;%单个锚链长度
m5=28.12*L5;%单个锚链质量
p5=7850;%钢的密度
V5=m5/p5;%单个锚链体积
G5=m5*g;
B5=p*g*V5;T5=zeros(nnn+1,1);
si5=zeros(nnn+1,1);
a5=zeros(nnn+1,1);
h5=zeros(nnn+1,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:nnnsi5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end
%% 绘图6
x=zeros(nnn+1,1);
y=zeros(nnn+1,1);
x(1)=0;
y(1)=0;
count=nnn+1;
count00=0;%链接漂浮个数
for ii=1:nnnif a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(nnn+1,1);%总高度
r0=r1+sum(r2)+r3+x(nnn+1,1);%游动半径
HH(count0)=H;
end
% plot(x,y,'linewidth',2)
% axis([0 20 0 20])
% xlabel('x','FontSize',28);
% ylabel('y','FontSize',28);
% set(gca,'FontSize',28,'linewidth',1);

3.3 代入遍历所得浮标吃水深度以及锚链型号与长度以及重物球的质量MATLAB程序

%% 浮标1
g=9.807;%重力加速度
p=1025;%海水密度
m1=1000;%浮标质量
v=12;%海面风速
v0=1.5;%海水水速
nnn=125;%锚链节个数h1=1.841;%吃水深度m1=1000;%浮标质量
V1=pi*h1;%浮标体积
B1=p*g*V1;%浮标最大浮力
G1=m1*g;%浮标重力
S1=2*h1;%受水投影面积
Ffeng=0.625*((2-h1)*2)*v^2;%F风
Fshui1=374*S1*v0^2;si1=atand((Fshui1+Ffeng)/(B1-G1));
T1=sqrt((B1-G1)^2+(Fshui1+Ffeng)^2);r1=0;%横坐标长度%% 钢管2
m2=10;%每节钢管质量
V2=pi*0.025^2*1;%钢管体积
L2=1;
G2=m2*g;
B2=p*g*V2;T2=zeros(5,1);
si2=zeros(5,1);
a2=zeros(4,1);
T2(1)=T1(1);si2(1)=si1(1);
S2=0.05*1;
Fshui2=374*S2*v0^2;
for n=1:4si2(n+1)=atand((T2(n)*sind(si2(n))+Fshui2)/(T2(n)*cosd(si2(n))+B2-G2));T2(n+1)=sqrt((T2(n)*cosd(si2(n))+B2-G2)^2+(T2(n)*sind(si2(n))+Fshui2)^2);a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
end
% for i=1:3
%     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
%     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
%     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
%     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
%     sol2=struct2cell(sol2);%结构体转元胞数组
%
%     index2=find(sol2{1,1}>0);%寻找大于0的解
%     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
%     si2(i+1)=double(sol2{2,1}(index2(1),1));
%     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
%     fT2(i+1)=T2(i+1);
%     fsi2(i+1)=si2(i+1);
%     fa2(i+1)=a2(i+1);
% end
h2=L2*cosd(a2);%高度
r2=abs(L2*sind(a2));%横坐标长度%% 重物球4
m4=5000;%重物球质量
p4=7850;%钢的密度
V4=m4/p4;%重物球体积
G4=m4*g;
B4=p*g*V4;
T41=G4-B4;%% 钢桶3
m3=100;%钢桶质量
V3=pi*0.15^2*1;%钢桶体积
L3=1;
G3=m3*g;
B3=p*g*V3;T31=T2(5);%取第2段的值
si31=si2(5);
S3=0.3*1;
Fshui3=374*S3*v0^2;
T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31)+Fshui3)^2);
si32=atand((T31*sind(si31)+Fshui3)/(T31*cosd(si31)+B3-G3-T41));syms fa3
eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
[fa3]=solve(eq3,fa3);%解方程
a3=double(fa3);
index3=find(a3>0);
a3=double(a3(index3));h3=L3*cosd(a3);%高度
r3=abs(L3*sind(a3));%横坐标长度%% 各锚链节5
L5=0.180;%单个锚链长度
m5=28.12*L5;%单个锚链质量
p5=7850;%钢的密度
V5=m5/p5;%单个锚链体积
G5=m5*g;
B5=p*g*V5;T5=zeros(nnn+1,1);
si5=zeros(nnn+1,1);
a5=zeros(nnn+1,1);
h5=zeros(nnn+1,1);
T5(1)=T32;si5(1)=si32;%取第3段的值
a5(1)=si32;
for i=1:nnnsi5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));if si5(i+1)<0abs(si5(i+1));endif abs(a5(i+1))<1a5(i+1)=90;end
end
%% 绘图6
x=zeros(nnn+1,1);
y=zeros(nnn+1,1);
x(1)=0;
y(1)=0;
count=nnn+1;
count00=0;%链接漂浮个数
for ii=1:nnnif a5(count)==0a5(count)=90;endx(ii+1)=x(ii)+L5*sind(a5(count));y(ii+1)=y(ii)+L5*cosd(a5(count));count=count-1;if a5(ii)~=90count00=count00+1;end
end
H=h1+sum(h2)+h3+y(nnn+1,1);%总高度
r0=r1+sum(r2)+r3+x(nnn+1,1);%游动半径
plot(x,y,'linewidth',2)
axis([0 20 0 20])
xlabel('x','FontSize',28);
ylabel('y','FontSize',28);
set(gca,'FontSize',28,'linewidth',1);
% %% 浮标1
% g=9.807;%重力加速度
% p=1025;%海水密度
% m1=1000;%浮标质量
% v=24;%海面风速
%
% h1=0.697;%吃水深度
%
% V1=pi*h1;%浮标吃水体积
% syms fT1;%拉力T（下同）
% syms fsi1;%角度SITA（下同）
%
% G1=m1*g;
% B1=p*g*V1;
% Ffeng=0.625*((2-h1)*2)*v^2;
% eq11=B1-G1-fT1*cosd(fsi1);
% eq12=Ffeng-fT1*sind(fsi1);
% [fT1,fsi1]=solve(eq11,eq12,fT1,fsi1);
% index1=find(fT1>0);
% T1=double(fT1(index1));
% si1=double(fsi1(index1));
% r1=0;%横坐标长度
%
% %% 钢管2
% m2=10;%每节钢管质量
% V2=pi*0.025^2*1;%钢管体积
% L2=1;
% G2=m2*g;
% B2=p*g*V2;
%
% T2=zeros(5,1);
% si2=zeros(5,1);
% a2=zeros(4,1);
% T2(1)=T1(1);si2(1)=si1(1);
% for n=1:4
%     si2(n+1)=atand((T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+B2-G2));
%     T2(n+1)=(T2(n)*sind(si2(n)))/(sind(si2(n+1)));
%     a2(n)=atand((T2(n+1)*sind(si2(n+1))+T2(n)*sind(si2(n)))/(T2(n)*cosd(si2(n))+T2(n+1)*cosd(si2(n+1))));
% end
% % for i=1:3
% %     eq21=fT2(i)*cosd(fsi2(i))+B2-m2*g-fT2(i+1)*cosd(fsi2(i+1));%受力
% %     eq22=fT2(i)*sind(fsi2(i))-fT2(i+1)*sind(fsi2(i+1));%受力
% %     eq23=0.5*fT2(i)*L2*sind(fa2(i+1)-fsi2(i))-0.5*fT2(i+1)*L2*sind(fsi2(i+1)-fa2(i+1));%力矩
% %     sol2=solve(eq21,eq22,eq23,fT2(i+1),fsi2(i+1),fa2(i+1));%解方程
% %     sol2=struct2cell(sol2);%结构体转元胞数组
% %
% %     index2=find(sol2{1,1}>0);%寻找大于0的解
% %     T2(i+1)=double(sol2{1,1}(index2(1),1));%写入答案
% %     si2(i+1)=double(sol2{2,1}(index2(1),1));
% %     a2(i+1)=double(sol2{3,1}(index2(1),1))+180;
% %     fT2(i+1)=T2(i+1);
% %     fsi2(i+1)=si2(i+1);
% %     fa2(i+1)=a2(i+1);
% % end
% h2=L2*cosd(a2);%高度
% r2=abs(L2*sind(a2));%横坐标长度
%
% %% 重物球4
% m4=1200;%重物球质量
% p4=7850;%钢的密度
% V4=m4/p4;%重物球体积
% G4=m4*g;
% B4=p*g*V4;
% T41=G4-B4;
%
% %% 钢桶3
% m3=100;%钢桶质量
% V3=pi*0.15^2*1;%钢桶体积
% L3=1;
% G3=m3*g;
% B3=p*g*V3;
%
% T31=T2(5);%取第2段的值
% si31=si2(5);
% T32=sqrt((T31*cosd(si31)+B3-G3-T41)^2+(T31*sind(si31))^2);
% si32=atand((T31*sind(si31))/(T31*cosd(si31)+B3-G3-T41));
% syms fa3
% eq3=0.5*T32*L3*sind(si32-fa3)-0.5*T41*L3*sind(fa3)-0.5*T31*L3*sind(fa3-si31);%力矩平衡
% [fa3]=solve(eq3,fa3);%解方程
% a3=double(fa3);
% index3=find(a3>0);
% a3=double(a3(index3));
%
% h3=L3*cosd(a3);%高度
% r3=abs(L3*sind(a3));%横坐标长度
%
% %% 各锚链节5
% L5=0.105;%单个锚链长度
% m5=7*L5;%单个锚链质量
% p5=7850;%钢的密度
% V5=m5/p5;%单个锚链体积
% G5=m5*g;
% B5=p*g*V5;
%
% T5=zeros(nnn+1,1);
% si5=zeros(nnn+1,1);
% a5=zeros(nnn+1,1);
% h5=zeros(nnn+1,1);
% T5(1)=T32;si5(1)=si32;%取第3段的值
% a5(1)=si32;
% for i=1:nnn
%     si5(i+1)=atand((T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+B5-G5));
%     T5(i+1)=sqrt((T5(i)*cosd(si5(i))+B5-G5)^2+(T5(i)*sind(si5(i)))^2);
%     a5(i+1)=atand((T5(i+1)*sind(si5(i+1))+T5(i)*sind(si5(i)))/(T5(i)*cosd(si5(i))+T5(i+1)*cosd(si5(i+1))));
%     if si5(i+1)<0
%         abs(si5(i+1));
%     end
%     if abs(a5(i+1))<1
%        a5(i+1)=90;
%     end
% end
% %% 绘图6
% x=zeros(nnn+1,1);
% y=zeros(nnn+1,1);
% x(1)=0;
% y(1)=0;
% count=nnn+1;
% count00=0;%链接漂浮个数
% for ii=1:nnn
%     if a5(count)==0
%        a5(count)=90;
%     end
%     x(ii+1)=x(ii)+L5*sind(a5(count));
%     y(ii+1)=y(ii)+L5*cosd(a5(count));
%     count=count-1;
%     if a5(ii)~=90
%         count00=count00+1;
%     end
% end
% H=h1+sum(h2)+h3+y(nnn+1,1);%总高度
% r0=r1+sum(r2)+r3+x(nnn+1,1);%游动半径
% hold on
% plot(x,y,'linewidth',2)
% axis([0 20 0 20])
% xlabel('x','FontSize',28);
% ylabel('y','FontSize',28);
% set(gca,'FontSize',28,'linewidth',1);
%
% legend('悬链线模型结果','受力递推模型结果');
%
%% 导出数据7
filterName=["钢管角度","钢桶角度","浮标吃水深度","浮标游动半径","锚链末端夹角","链结漂浮个数","系统总高度","重物质量","锚链型号"];
xlswrite('problem3.xlsx',filterName,'v=12','A1');
xlswrite('problem3.xlsx',a2,'v=12','A2');
xlswrite('problem3.xlsx',a3,'v=12','B2');
xlswrite('problem3.xlsx',h1,'v=12','C2');
xlswrite('problem3.xlsx',r0,'v=12','D2');
xlswrite('problem3.xlsx',90-a5(nnn),'v=12','E2');
xlswrite('problem3.xlsx',count00,'v=12','F2');
xlswrite('problem3.xlsx',H,'v=12','G2');
xlswrite('problem3.xlsx',m4,'v=12','H2');
xlswrite('problem3.xlsx',5,'v=12','I2');

全国大学生数学建模竞赛2016A题系泊系统的设计MATLAB程序相关推荐

  1. 2021 年高教社杯全国大学生数学建模竞赛A题分析

    2021 年高教社杯全国大学生数学建模竞赛A题分析 题目 赛题分析 前言 问题一分析 问题二分析 问题三分析 题目 A 题 "FAST"主动反射面的形状调节 中国天眼--500 米 ...

  2. 2020年高教社杯全国大学生数学建模竞赛 C题思路

    2020年高教社杯全国大学生数学建模竞赛 C题 中小微企业的信贷决策 本文旨在为广大热爱建模的朋友们提供2020年数学建模C题的思路和解法. 问题回顾 在实际中,由于中小微企业规模相对较小,也缺少抵押 ...

  3. 2020年全国大学生数学建模竞赛B题穿越沙漠问题——建立整数线性规划模型（ILP）——通过LINGO求解

    2020年全国大学生数学建模竞赛B题 穿越沙漠 题目是讲玩家在不同地图下穿越沙漠,所获得的资金数要最多(大概是这个意思).然后通过文章的描述又总结了N个约束条件.整体的思路就是对资金最大化作为目标函数 ...

  4. 2020年高教社杯全国大学生数学建模竞赛C题 第一问详细解答+代码

    2020年高教社杯全国大学生数学建模竞赛C题 第一问详细解答+代码 本文摘自小编自己的参赛论文与经历,小编获得了2020年高教社杯国奖,有问题的同学们可私聊博主哦. 1. 问题分析 问题一主要围绕信贷 ...

  5. 【数学建模】2003年全国大学生数学建模竞赛B题求解

    目录 [数学建模]2003年全国大学生数学建模竞赛B题求解 [数学建模]2003年全国大学生数学建模竞赛B题求解 model: title CUMCM-2003B-01; sets: cai / 1. ...

  6. 2022年全国大学生数学建模竞赛C题思路

    一.思路获取方式 获取代码方式: 2022年全国大学生数学建模竞赛赛题思路 备注: 点击上面蓝色字体2022年全国大学生数学建模竞赛赛题思路,扫描上面二维码,付费29.9元订阅海神之光博客付费专栏20 ...

  7. 2022年全国大学生数学建模竞赛E题思路

    一.思路获取方式 获取代码方式: 2022年全国大学生数学建模竞赛赛题思路 备注: 点击上面蓝色字体2022年全国大学生数学建模竞赛赛题思路,扫描上面二维码,付费29.9元订阅海神之光博客付费专栏20 ...

  8. 2022年全国大学生数学建模竞赛A题思路

    一.思路获取方式 获取代码方式: 2022年全国大学生数学建模竞赛赛题思路 备注: 点击上面蓝色字体2022年全国大学生数学建模竞赛赛题思路,扫描上面二维码,付费29.9元订阅海神之光博客付费专栏20 ...

  9. 1998年全国大学生数学建模竞赛A题——投资的收益和风险数模P133|lingo，matlab

    1998年全国大学生数学建模竞赛A题 目录 题目 问一 用lingo求解 用matlab求解 问2 题目 市场上有n 种资产(如股票.债券.-)Si ( i=1,-n) 供投资者选择,某公司有数额为M ...

  10. 尖峰法聚类：2021 年高教社杯全国大学生数学建模竞赛 E题 中药材的鉴别 问题1

    PeakCluster是Lu优化库中的一个函数,该函数利用数据曲线尖峰形状和位置进行聚类分析. 例子:2021 年高教社杯全国大学生数学建模竞赛 E题 中药材的鉴别 问题1:根据附件 1 中几种药材的 ...

最新文章

  1. pytorch 测试每一类_DeepFM全方面解析（附pytorch源码）
  2. js_xpath_搞不定的东西
  3. UDP数据转发解决WiFi与有限以太网之间控制命令传递：RGBLink
  4. display: none;、visibility: hidden、opacity=0区别总结
  5. Jacobi迭代法与Gauss-Seidel迭代法
  6. oracle中lock和latch的用途
  7. Bitmap的内存占用
  8. Aspx 页面生命周期
  9. 【文末有福利】吸烟致癌，是基因的错吗？
  10. CV Code|计算机视觉开源周报20200602期～文末送书
  11. ImageMagick还是GraphicsMagick?
  12. docker compose详解
  13. 6003.mavlink协议自定义消息编程
  14. java零碎要点012---linux Centos下编译、运行、调试java程序
  15. 应用程序框架实战十八:DDD分层架构之聚合
  16. 中英文对照 —— 软件与病毒、电子与硬件
  17. 创建系统镜像_学会这些, 操作docker image镜像就够了！
  18. 八大排序算法-java实现
  19. java 无驱动socket连接热敏小票打印机示例，编写自定义模板 芯烨/xprinter，附工具类即开即用
  20. esxi虚拟化服务器端口聚合,配置链路聚合组处理分布式端口组的流量

热门文章

  1. python写入access_使用Python对Access读写操作方法详解
  2. 计算机财务管理中的函数,浅析几种常用的Excel函数在财务管理中的应用
  3. 淡雅简洁商业汇演商业计划书PPT模板
  4. 防止链接和二维码被微信拦截（被封锁、被屏蔽、被和谐）的最新方法——MaxJump
  5. ## python爬取MM131整站图片到本地
  6. 【数学建模】论文模板和latex模板
  7. 火力全开，同时分解(切脸)多个视频
  8. android+p开机动画,android开机动画[转]
  9. 海康威视4路播放封装----安卓开发
  10. 庖丁解D,游刃有余---Discuz!免费版安全性分析(转)