2021高教杯数学建模A

数学建模比赛A题第二问程序分享

clear all;clc;close all;

%载入节点数据
% [NodeLoc,NodeId1]=xlsread('附件1.csv');
% [NodeLoc2,NodeId2]=xlsread('附件2.csv');
% [~,NodeId3]=xlsread('附件3.csv');
% save NodeInfo NodeLoc NodeId1 NodeLoc2 NodeId2 NodeId3
load NodeInfo


NodeNum=size(NodeLoc,1);     %NodeNum=NodeLoc矩阵第一列的长度
R=abs(min(NodeLoc(:,3)));    %R=NodeLoc矩阵第三列的最小值的绝对值
f=0.466;                     %定义常量f=0.466
MaxMove=0.6;

%对坐标系进行旋转
rad=pi/180;          %rad=将一弧度分成180份（1°对应的弧度）
alpha= 36.795*rad;   %alpha=36.795°所对应的弧度（将方位角化为弧度）
beta= (90-78.169)*rad;  %beta=仰角的补角且将该角弧度化
rat_mat1=[cos(alpha) sin(alpha) 0 0;-sin(alpha) cos(alpha) 0 0;0 0 1 0;0 0 0 1;];
%写成矩阵的形式更加直观一些
%将坐标轴绕z轴旋转alpha弧度
%        cos(alpha)   sin(alpha)   0   0
%       -sin(alpha)   cos(alpha)   0   0
%            0            0        1   0
%            0            0        0   1
rat_mat2=[cos(beta) 0 -sin(beta) 0;0 1 0 0;sin(beta)  0 cos(beta) 0 ;0 0 0 1;];
%将坐标绕y轴旋转beta弧度
%        cos(beta)   0   -sin(beta)   0
%            0       1       0        0
%        sin(beta)   0    cos(beta)   0
%            0       0       0        1
rat_mat2I=[cos(-beta) 0 -sin(-beta) 0;0 1 0 0;sin(-beta)  0 cos(-beta) 0 ;0 0 0 1;];
%将坐标轴绕y轴旋转-beta弧度
rat_mat1I=[cos(-alpha) sin(-alpha) 0 0;-sin(-alpha) cos(-alpha) 0 0;0 0 1 0;0 0 0 1;];
%将坐标轴绕z轴旋转-alpha弧度
ratmat12=rat_mat1*rat_mat2;     %将坐标轴先绕z轴旋转alpha再绕y轴旋转beta
ratmat12I=rat_mat2I*rat_mat1I;  %将坐标轴先绕y轴旋转-beta再绕z轴旋转-alpha
NodeLocr=NodeLoc;NodeLocr(:,3)=NodeLocr(:,3)+R;NodeLocr(:,4)=1;  %先将NodeLoc赋值给NodeLocr,然后将NodeLocr第三列各加R，再将该矩阵加一列全是1的第四列
NodeLocr=NodeLocr*ratmat12;     %将NodeLocr这个矩阵经过ratmat12旋转（将原来的坐标加了R又加了一列1后再进行了旋转）

[~,Indmid]=min(NodeLocr(:,3));  %找到NodeLocr矩阵中第三列的最小值然后赋值给Indmid
NodeLocOffset=NodeLocr(Indmid,1:3);   %找到NodeLocr矩阵的第Indmid行然后赋值给NodeLocr0ffset
NodeLocr(:,1)=NodeLocr(:,1)-NodeLocOffset(1);  %将NodeLocr的第一列都减去NodeLocr0ffset的第一个值
NodeLocr(:,2)=NodeLocr(:,2)-NodeLocOffset(2);  %将NodeLocr的第二列都减去NodeLocr0ffset的第二个值
NodeLocr(:,3)=NodeLocr(:,3)-NodeLocOffset(3)-R;%将NodeLocr的第三列都减去NodeLocr0ffset的第三个值再减R
%总结一下：以上这部分坐标转换是先将原来节点的坐标经过坐标旋转后，再找出旋转后的坐标最低点（如果画图的话就是旋转后图像Z最小的那个坐标）
%          然后将原来的每个点减去这个最低点的坐标（相当于是将原点交给了最低点），这样就实现了平移的效果，到这里坐标的变换就结束了

MeshNum=size(NodeId3,1)-1;   %提取NodeId3矩阵第一列的长度-1(看看NodrId3矩阵的值就知道为什么要-1了）赋值给变量MeshNum
MeshGrid=zeros(MeshNum,3);   %生成一个MeshNum乘3阶的全零矩阵MeshGrid
for ii=1:NodeNum
    NodeId1Name=NodeId1{ii+1,1};
    Row=NodeId1Name(1)-'A'+1;
    Col=str2double(NodeId1Name(2:end))+1;
    Node_ID(Row,Col)=ii;
end
%这个循环结构的作用是将附件一的第一列分出来，然后将每个元素的字母和数字再分开，将字母转化为数字然后作为行，数字作为列，将第一列
%前的序号排进矩阵

for ii=1:MeshNum
    for jj=1:3
        NodeId1Name=NodeId3{ii+1,jj};
        Row=NodeId1Name(1)-'A'+1;
        Col=str2double(NodeId1Name(2:end))+1;
        MeshGrid(ii,jj)=Node_ID(Row,Col);
    end
end
 %这个循环结构的作用是将附件三以上一个循环的形式将序号保存到矩阵中
 
TargetLoc=NodeLocr(:,1:3); 
RealLoc=NodeLocr(:,1:3);
moveDisR=zeros(1,NodeNum);   %生成一个1乘NodeNum的全零矩阵
for ii=1:NodeNum
    LocOld=NodeLocr(ii,1:3);  %从转换后的坐标中拿出一组坐标
    [LocNew]=locSfun(LocOld,R,f);   %通过自定义函数计算出相对应的理想抛物面上那个节点的坐标
    TargetLoc(ii,:)=LocNew;
    moveDis=norm(LocOld-LocNew);   %计算基态节点与理想抛物面上相对应节点之间的距离
    if (moveDis>0)
        MovDir=LocNew-LocOld;
        MovDir=MovDir/norm(MovDir);  %计算出基态节点和抛物面对应节点连线对于坐标轴的余弦值
        
        if (moveDis<MaxMove)
            RealLoc(ii,:)=LocNew;    %如果距离小于可移动范围，则移动的距离就是相对应节点的距离
        else
            moveDis=0.6;
            RealLoc(ii,:)=LocOld+MovDir*MaxMove;  %如果大于或者等于最大移动距离，则将将最大移动距离分解到节点的连线方向再加上基态坐标
        end
        %这个循环结构的作用是计算出了考虑移动条件下的实际接受抛物面上各节点相对应的坐标
        
        LocDir=LocOld/norm(LocOld);
        movesign=-sign(MovDir*LocDir');
        moveDisR(ii)=movesign*moveDis;   %计算实际促动器的伸缩量
    end
end
%//
TargetLocr=TargetLoc;
TargetLocr(:,3)=TargetLocr(:,3)+R+NodeLocOffset(3);
TargetLocr(:,2)=TargetLocr(:,2)+NodeLocOffset(2);
TargetLocr(:,1)=TargetLocr(:,1)+NodeLocOffset(1);
TargetLocr(:,4)=1;TargetLocr=TargetLocr*ratmat12I;
TargetLocr(:,3)=TargetLocr(:,3)-R;
RealLocr=RealLoc;
RealLocr(:,3)=RealLocr(:,3)+R+NodeLocOffset(3);
RealLocr(:,2)=RealLocr(:,2)+NodeLocOffset(2);
RealLocr(:,1)=RealLocr(:,1)+NodeLocOffset(1);
RealLocr(:,4)=1;RealLocr=RealLocr*ratmat12I;
RealLocr(:,3)=RealLocr(:,3)-R;
%//这一部分是进行逆向旋转坐标轴后的结果
 
%绘图调色
figure(1);clf;hold on;
plot3(NodeLoc(:,1),NodeLoc(:,2),NodeLoc(:,3),'r.');
plot3(TargetLocr(:,1),TargetLocr(:,2),TargetLocr(:,3),'k.');
plot3(RealLocr(:,1),RealLocr(:,2),RealLocr(:,3),'b.');
legend('基准态','理想抛物面','调整后');
axis([-320 320 -320 320 -320 0]);
axis equal

figure(2);clf;hold on;
PlotModel(NodeLoc,MeshGrid,'b');
PlotModel(TargetLocr,MeshGrid,'r');
title('基准态（蓝色）与理想抛物面（红色）')

figure(3);clf;hold on;
PlotModel(NodeLoc,MeshGrid,'b');
PlotModel(RealLocr,MeshGrid,'r');
title('基准态（蓝色）与调整后抛物面（红色）')

figure(4);clf;hold on;
PlotModel(TargetLocr,MeshGrid,'b');
PlotModel(RealLocr,MeshGrid,'r');
title('理想抛物面（蓝色）与调整后抛物面（红色）');

%导出附件4
fid=fopen('附件4.csv','wt');
fprintf(fid,'对应主索节点编号,伸缩量（米）\n');
for ii=1:NodeNum
   fprintf(fid,'%s,',NodeId1{ii+1,1});
   fprintf(fid,'%0.4f,\n',moveDisR(ii));  
end
fclose(fid);
save res2