数学建模2023-A太阳镜厂代码认识

  P=np.column_stack((p1,new_column))  #得到每个镜子的x,y,z序列

?

nl=nl/np.linalg.norm(nl)  #得到单位法向量

?

 for dx in np.arange(-W/2,W/2+0.1,delta_t):

?

 indices_in_circle=np.where(Dis[:,i]==1)[0] #取周围半径

?

 Di_b=Tb.T.dot(Di_d-B) #A镜上的点 从地面坐标系->B镜坐标系

?

总代码

1、

import numpy as np
import pandas as pd
from math import cos,sin,acos,asin,pi,exp,sqrt
import time
P=pd.read_excel(r'附件.xlsx').values
Dis=np.load('distance.npy')

ST_0=[9,10.5,12,13.5,15]
D_0=[337,0,31,61,92]
delta_t=1.5
x_c=0;y_c=0;L=W=6
H0=80;H1=8;H=4
HR=3.5

def F(L,W,H,D_0,ST_0,delta_t,H0,H1,HR,Dis,p1,x_c,y_c):
    N=(W/delta_t+1)**2*4
    yita_ref=0.92  #镜面反射率
    S=L*W

    new_column=np.array([H for i in range(len(p1))]) #安装高度4m
    P=np.column_stack((p1,new_column))  #得到每个镜子的x,y,z序列
    #################初始化yita########################
    Yita=np.zeros([len(ST_0),len(P)])
    Yita_trunc=np.zeros([len(ST_0),len(P)])
    Yita_at=np.zeros([len(ST_0),len(P)])
    Yita_cos =np.zeros([len(ST_0),len(P)])
    Yita_sb =np.zeros([len(ST_0),len(P)])

    E_0=np.zeros([len(D_0),len(ST_0)])
    Y1=Y2=Y3=Y4=np.zeros([len(D_0),len(ST_0)])

    DF=np.zeros(len(D_0),5)

    for di, D in enumerate(D_0):
        for sti,ST in enumerate(ST_0):
           #计算太阳位置 以及相关参数
           fai=39.4*pi/180   #得到用弧度表示的当地纬度
           delta=asin(sin(2*pi*D/365)*sin(2*pi*23.45/360))   #delta为太阳赤纬角
           w=(pi/12)*(ST-12)   #得到为太阳时角

           sin_as=cos(delta)*cos(fai)*cos(w)+sin(delta)*sin(fai)#太阳高度角
           cos_as=sqrt(1-sin_as**2)

           cos_rs=(sin(delta)-sin_as*sin(fai))/(sqrt(1-sin_as**2)*cos(fai)) #太阳方位角
           if abs(cos_rs)>1:     #一般是不会出现>1的情况，若出现了，说明值有一些问题
               cos_rs=cos_rs/abs(cos_rs)

           if ST<=12:
               sin_rs=sqrt(1-cos_rs**2)
           else:
               sin_rs=-sqrt(1-cos_rs**2)

           a=0.4237-0.00821-(6-3)**2
           b=0.5055+0.00595-(6.5-3)**2
           c=0.2711+0.01858*(2.5-3)**2
           DNI=1.366*(a+b*exp(-c/sin_as))    #得到法向直接辐射辐照度 DNI
           A0=np.array([x_c,y_c,H0])#集热器中心
           Ls=np.array([-cos_as*sin_rs,-cos_as*cos_rs,-sin_as]) #入射光线的方向向量
           s=np.array([1,0,0])

           for i in range(len(P)):
               #A镜中点  反射向量
               Di=P[i]
               Lr=A0-Di
               nl=-Ls+Lr/np.linalg.norm(Lr) #A的法向量
               nl=nl/np.linalg.norm(nl)  #得到单位法向量

               beta1=asin(nl.dot([0,0,1])/np.linalg.norm(nl)) #俯仰角

               n0=np.array([nl[0],nl[1],0])
               if nl[1]>=0:
                   alpha1=acos(n0.dot(s)/np.linalg.norm(n0)) #方位角
               else:
                   alpha1 = -acos(n0.dot(s) / np.linalg.norm(n0))

               #A的旋转矩阵
               Ta=np.array([
                   [cos(alpha1)*cos(pi/2-beta1),-sin(alpha1),cos(alpha1)*sin(pi/2-beta1)],
                   [sin(alpha1) * cos(pi / 2 - beta1), cos(alpha1), sin(alpha1) * sin(pi / 2 - beta1)],
                   [-sin(pi/2-beta1),0,cos(pi/2-beta1)]
               ])

               light=0
               empty=0
               barr_tower=0
               barr_s=0
               barr_r=0
               #遍历每一个点
               for dx in np.arange(-W/2,W/2+0.1,delta_t):
                   for dy in np.arange(-L/2,L/2+0.1,delta_t):
                       Dxy=np.array([dx,dy,0]) #A镜上的某个点在A镜坐标系
                       Di_d=Ta.dot(Dxy)+Di #A镜上的点 转置到地面坐标系
                       ##############遍历每一个入射光线圆锥的光线##############
                       #for the1 in np.arange(0.001,0.00465,0.002):
                       the1=0.02
                       for the2 in np.arange(0,2*pi,pi/2):
                            if_barr=0
                            #g是入射光线在主光线锥体系的坐标
                            g=np.array([sin(the1)*cos(the2),
                                      sin(the1)*sin(the2),
                                      cos(the1)])
                            #根据入射主光线，计算主光线锥体系的旋转矩阵 Ls：入射的太阳光线
                            v=pi/2-acos(Ls.dot(np.array([0,0,1]))/np.linalg.norm(Ls))  #锥体光线的俯仰角
                            nl_g0=np.array([Ls[0],Ls[1],0])
                            if Ls[1]>=0:
                                   u=acos(nl_g0.dot(s)/np.linalg.norm(nl_g0)) #锥体主光线方位角
                            else:
                                   u = -acos(nl_g0.dot(s) / np.linalg.norm(nl_g0)) #锥体主光线方位角
                            T_s=np.array([
                                    [cos(u)*cos(pi/2-v),-sin(u),cos(u)*sin(pi/2-v)],
                                    [sin(u) * cos(pi / 2 - v), cos(u), sin(u) * sin(pi / 2 - v)],
                                    [-sin(pi/2-v),0,cos(pi/2-v)]
                                    ])   #光锥到大地的转换矩阵
                            g_d=T_s.dot(g) #转置到地面坐标系 g:入射光锥中的某条入射线
                            g_r=g_d-2*g_d.dot(nl)*(nl) #对应的反射向量  #nl:A的法向量

                            ##################一、判断入射光线是否被遮挡##########################
                            a,b,c=g_d #入射光线的x、y、z
                            x0,y0,z0=Di_d  #镜子上的点的坐标
                            delta_tower=4*(a*(x0-x_c)+b*(y0-y_c))**2-4*(a**2+b**2)*((x0-x_c)**2+(y0-y_c)**2-HR**2)
                            if delta_tower>=0: #计算两个交点根
                                   t1=(-2*(a*(x0-x_c)+b*(y0-y_c))+sqrt(delta_tower))/(2*(a**2+b**2))
                                   t2 = (-2 * (a * (x0 - x_c) + b * (y0 - y_c)) - sqrt(delta_tower)) / (  2 * (a ** 2 + b ** 2))
                            if min(t1*c+z0,t2*z0)<=(H0+H1/2)and min(t1*c+z0,t2*z0)>=0:
                                   barr_tower+=1
                                   continue

                               #########二、判断入射光线是否被其他光镜遮挡##############
                            indices_in_circle=np.where(Dis[:,i]==1)[0] #取周围半径
                            for j in indices_in_circle: #len(P)
                                   if i==j:
                                       continue
                                   #B的中心坐标 alpha,beta直接调取提前计算的
                                   B=P[j]
                                   Lrb=A0-B  #反射光线
                                   nlb=-Ls+Lrb/np.linalg.norm(Lrb) #法向量
                                   #beta2,alpha2=all_alpha_beta[j,:]

                                   beta2=asin(nlb.dot[0,0,1])/np.linalg.norm(nlb) #俯仰角

                                   n0b=np.array(nlb[0],nlb[1],0)
                                   if nlb[1]>=0:
                                       alpha2=acos(n0b.dot(s)/np.linalg.norm(n0b) )#方位角
                                   else:
                                       alpha2 = -acos(n0b.dot(s) / np.linalg.norm(n0b))

                                   Tb=np.array([
                                    [cos(alpha2)*cos(pi/2-beta2),-sin(alpha2),cos(alpha2)*sin(pi/2-beta2)],
                                    [sin(alpha2) * cos(pi / 2 - beta2), cos(alpha2), sin(alpha2) * sin(pi / 2 - beta2)],
                                    [-sin(pi/2-beta2),0,cos(pi/2-beta2)]
                                    ])

                                   Di_b=Tb.T.dot(Di_d-B) #A镜上的点 从地面坐标系->B镜坐标
                                   g_b=Tb.T.dot(g_d) #A入射光线 从地面坐标系->B镜坐标系

                                   t=-Di_b[2]/g_b[2] #计算入射光线的在B镜坐标系的交点
                                   x_b=g_b[0]*t+Di_b[0]
                                   y_b = g_b[1] * t + Di_b[1]
                                   D0=np.array([x_b,y_b,0])
                                   D0=Tb.dot(D0)+B #交点转到地面

                                   if abs(x_b)<=W/2 and abs(y_b)<=L/2 and D0[2]>Di_d[2]: #如果被遮挡了，就直接算下一条入射的线
                                       barr_s+=1
                                       if_barr=1
                                       break

                                   ##################三、对应的反射###############
                                   g_r=g_d-2*g_d.dot(nl)*(nl) #对应的反射向量 #nl:A的法向量
                                   g_r_b=Tb.T.dot(g_r) #转置到b镜面坐标系

                                   t=-Di_b[2]/g_r_b[2] #计算反射光线的交点
                                   x_b=g_r_b[0]*t+Di_b[0]
                                   y_b = g_b[1] * t + Di_b[1]

                                   D0=np.array([x_b,y_b,0])
                                   D0 = Tb.dot(D0) + B  # 交点转到地面

                                   if abs(x_b)<=W/2 and abs(y_b)<=L/2 and D0[2]>Di_d[2]:
                                       barr_r+=1
                                       if_barr=1
                                       break

                            ###########四、是否吸收################
                            if if_barr==0:
                                   a,b,c=g_r
                                   x0,y0,z0=Di_d
                                   delta_recieve=4*(a*(x0-x_c)+b*(y0-y_c))**2-4*(a**2+b**2)*((x0-x_c)**2+(y0-y_c)**2-HR**2)
                                   if delta_recieve>=0:
                                       t1=( -2*(a*(x0-x_c)+b*(y0-y_c))+sqrt(delta_recieve))/(2*(a**2+b**2))
                                       t2=( -2*(a*(x0-x_c)+b*(y0-y_c))+sqrt(delta_recieve))/(2*(a**2+b**2))
                                       if   min(t1*c+z0,t2*c+z0)<=(H0+H1/2) and min(t1*c+z0,t2*c+z0)>=(H0-H1/2) :
                                           light+=1
                                       else:
                                           empty+=1
           #这里是这个时间点，计算第i个面板的数值，并记录，每个列表是1745的长度
           yita_sb=1-(barr_r+barr_s+barr_tower)/N
           yita_cos=abs(Ls.dot(-nl)/np.linalg.norm(Ls))
           HR0=np.linalg.norm(Ls)
           yita_at=0.99321-0.0001176*HR0+1.97e-8*(HR0**2)

           if N-barr_s-barr_r-barr_tower==0:
                 yita_trunc=1
           else:
                yita_trunc=(light)/(N-barr_s-barr_r-barr_tower)

           yita=yita_sb*yita_cos*yita_at*yita_trunc*yita_ref

           Yita_sb[sti,i]=yita_sb
           Yita_cos[sti, i] = yita_cos
           Yita_at[sti, i] = yita_at
           Yita_trunc[sti, i] = yita_trunc
           Yita[sti, i] = yita
        #这里是在D，ST的循环里，计算第sti个时间点
        E_0[di,sti]=DNI*sum(S*Yita[sti,:])

        Y1[di,sti]=np.mean(Yita[sti,:])
        Y2[di, sti] = np.mean(Yita_cos[sti, :])
        Y3[di, sti] = np.mean(Yita_sb[sti, :])
        Y4[di, sti] = np.mean(Yita_trunc[sti, :])
        #已经计算完一个具体时间点
        DF[di,0]=np.mean(Y1[di,:])
        DF[di, 1] = np.mean(2[di, :])
        DF[di, 2] = np.mean(Y3[di, :])
        DF[di, 3] = np.mean(Y4[di, :])
        DF[di, 4] = sum(E_0[di,:])/(len(P)*S)/len(ST_0)

        #print(DF)

    return DF
start_time=time.time()

Wp=F(L,W,H,D_0,ST_0,delta_t,H0,H1,HR,Dis,P,x_c,y_c)
print(Wp)

#记录程序运行时间
end_time=time.time()
elapsed_time_seconds=end_time-start_time
elapsed_time_minutes=round(elapsed_time_seconds/60,2)
print(f"代码运行时间：{elapsed_time_minutes}分钟")

2、

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import random
from math import cos,sin,acos,asin,pi,exp,sqrt,floor
import time
from Function import F

#导入题目的附件，并附加z值
p1=pd.read_excel(r's',header=None).values
HR=3.5
H=4
L=W=6
H0=80
H1=8

#导入每块板 在某时刻的 俯仰角与方位角
D_0=[275]
ST_0=[12]

#旋转角度
for gama in np.linspace(0,pi/3+0.01,4):
    TT=np.array([
        [cos(gama),-sin(gama)],
        [sin(gama),cos(gama)]
    ])
    p1=TT.dot(p1.T).T
    p1=p1

    Dis0=np.load('D')

    xy_len=range(-250,251,50)
    delta_t=1
    EP=np.zeros([len(xy_len),len(xy_len)])
    for m,x_c in enumerate(xy_len):
        for n,y_c in enumerate(xy_len):
            if x_c**2+y_c**2>250**2:
                continue
            EP0=0
            Pi=[]
            for i in range(len(p1)):
                dis=sqrt((p1[i,0]-x_c)**2+(p1[i,1]-y_c)**2)
                if dis>=100:
                    Pi.append(i)
            p1=p1[Pi]

            #生成新的连接矩阵
            Dis1=np.ones([len(Pi),len(Pi)])
            for i in range(len(Pi)):
                for j in range(len(Pi)):
                    Dis1[i,j]=Dis0[Pi[i],Pi[j]]

            m_n=floor(0.1*len(p1))
            for k in range(1,3):
                random_numbers=random.sample(list(range(len(p1))),m_n)
                p11=p1[random_numbers]
                Dis=np.ones([m_n,m_n])
                for i in range(m_n):
                    for j in range(m_n):
                        Dis[i,j]=Dis0[random_numbers[i],random_numbers[j]]

                Wp=F(L,W,H,D_0,ST_0,delta_t,H0,H1,HR,Dis,p11,x_c,y_c)
                EP0+=np.mean(Wp)*W*L*len(p1)/1000
            print(x_c,y_c,EP0/3)
            EP[m,n]=EP0

这段代码的作用是为了计算关于附件给定数据和某些参数的最优解。

具体来说，这份代码先在坐标系中旋转给定的数据，得到旋转后的坐标系，然后定义了一些常量和变量，包括板子的物理尺寸，初始的俯仰角和方位角，旋转角度等等。代码中还调用了自己编写的Function模块，这个模块里有一些函数，可以计算Wp和F，分别表示某段时间内WE算法的可用性和目标函数值。紧接着，代码中生成了一些新的矩阵并对其进行操作，最后输出计算出来的结果。

具体来说，循环中枚举了不同的旋转角度，通过旋转将数据映射到不同的坐标系中，然后在空间中枚举每个点(x_c, y_c)，其中点的坐标范围[-250, 250]。如果当前点与天空圆柱的交点不能连接到所有的板上，则无法继续计算此点的期望性能，直接跳过。接下来，计算可以到达所有板子的交点，这些点仅限于数据中离当前点距离小于100的点。这意味着我们不需要考虑从大距离发送的信号，因为这些信号方向不同，不能同时到达所有的板子。然后生成新的连接矩阵Dis1，并从中随机取出一定数量的点random_numbers，这些点的数量取决于数据中包含多少点。将选中的点的数据放入F函数中进行计算，计算结果存入EP矩阵中。最后输出EP矩阵，即为最终的计算结果。

?

3、

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import random
from math import cos,sin,acos,asin,pi,exp,sqrt,floor
import time
from Function_2 import F

#导入题目的附件，并附加z值
p1=pd.read_excel(r's',header=None).values
Dis0=np.load('s')

HR=3.5
H=4
L=W=6
H0=80
H1=8
delta_t=1.5 #步长

#导入每块板 在某时刻的 俯仰角与方位角
ST_0=[9,10.5,12,13.5,15]
D_0=[306]
time.time()

x_c=0
y_c=250
Pi=[]
start_time=time.time()

#判断r=100
for i in range(len(p1)):
    dis = sqrt((p1[i, 0] - x_c) ** 2 + (p1[i, 1] - y_c) ** 2)
    if dis >= 100:
        Pi.append(i)
p1 = p1[Pi]

Dis1=np.ones([len(Pi),len(Pi)])
for i in range(len(Pi)):
    for j in range(len(Pi)):
         Dis1[i,j]=Dis0[Pi[i],Pi[j]]

m_n=floor(0.3*len(p1))
random_numbers = random.sample(list(range(len(p1))), m_n)
p11 = p1[random_numbers]
Dis = np.ones([m_n, m_n])
for i in range(m_n):
    for j in range(m_n):
        Dis[i, j] = Dis0[random_numbers[i], random_numbers[j]]

Wp = F(L, W, H, D_0, ST_0, delta_t, H0, H1, HR, Dis, p11, x_c, y_c)

print(Wp)
#记录程序运行时间
end_time=time.time()
elapsed_time_seconds=end_time-start_time
elapsed_time_minutes=round(elapsed_time_seconds/60,2)
print(f"代码运行时间:{elapsed_time_minutes}分钟")

这段代码主要是在给定x_c和y_c的情况下，计算WE问题中WE算法的可用性Wp的值。

具体来说，代码使用了numpy和pandas等库加载了题目给定的附件，将其存储为数组p1和矩阵Dis0，并定义了一些常量和变量，如板子的物理尺寸L、W、H，以及时间间隔delta_t等。此外，代码还调用了自己编写的Function_2模块，并导入其中的F函数，F函数用来计算WE算法的可用性Wp。

接下来，代码对于每个点进行判断，判断距离该点最远的板子与该点之间的距离是否小于100，如果小于100，则记录该点的信息，并筛选掉距离该点最远的板子与该点之间距离大于等于100的点。然后，将筛选后的点生成一个新的矩阵Dis1，用于计算WE算法的可用性Wp。

再接下来，代码随机选取一定数目的点，并将这些点的数据放入F函数中进行计算，计算结果存入Wp中。最终输出Wp，即为计算结果。

最后代码输出了程序运行的时间。

?

4、

clc;clear;

k=2;   %计数
R=360
r=11;  %宽+5
n=fix(2*R/r);  %一排的点数

%第一排第一个点
W=[-R,R];
W1=W;
%计算第一排
for i =1:n
 W(k,1)=W(k-1,1)+r;
 W(k,2)=W(1,2);
 k=k+1;
end

%第二排第一个点
W(k,1)=-R+r/2;
W(k,2)=R-r*sqrt(3)/2;
W2=W(k,:);
k=k+1;

%计算第二排
for i =1:n
 W(k,1)=W(k-1,1)+r;
 W(k,2)=W(k-1,2);
 k=k+1;
end

n1=fix(2*R/(r*sqrt(3)/2))+1;
%从第三排开始
for i=3:n1
 if mod(i,2)==1 %计数行
   %先写该行的第一个点
   W(k,1)=W1(1,1);
   W(k,2)=W1(1,2)-sqrt(3)*r*(i-1)/2;
   k=k+1;

   for j=1:n
     W(k,1)=W(k-1,1)+r;
     W(k,2)=W(k-1,2);
     k=k+1;
   end
 else %偶数行
   %先写该行的第一个点
   W(k,1)=W2(1,1);
   W(k,2)=W2(1,2)-sqrt(3)*r*(i-2)/2;
   k=k+1;

   for j=1:n
     W(k,1)=W(k-1,1)+r;
     W(k,2)=W(k-1,2);
     k=k+1;
   end
 end
end

%计算范围内点
a=0;
for i=1:length(W)
 if W(i,1)^2+W(i,2)^2<=(350-(r-5)/2)^2%& W(i,1)^2+W(i,2)^2>=100^2
   a=a+1;
   W_e(a,:)=W(i,:);
 end
end

scatter(W_e(:,1),W_e(:,2),'.')

T=[cos(pi/4),-sin(pi/4);
  sin(pi/4),cos(pi/4)
];

WW=(T*W')';

%计算点距离
for i=1:length(W)
 for j=i:length(W)
   D(i,j)=sqrt((WW(i,1)-WW(j,1))^2+(WW(i,2)-WW(j,2))^2);
   D(j,i)=D(i,j);
 end
end

?

这段代码是用来生成WE问题中的点的位置。代码首先定义了一些参数和变量，如半径R、宽度r等。

然后，通过循环计算出了在圆形区域内的所有点的坐标位置，并将其存储在W中。

接下来，代码根据该圆的特殊特点，将得到的点进行旋转，旋转后存储在WW中。

最后，代码计算出点之间的距离矩阵D，用于后续计算WE算法中的Dis0。

最后，代码将生成的所有点进行可视化。

5、

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import random
from math import cos,sin,acos,asin,pi,exp,sqrt,floor
import time
from Function_3 import F

p1=pd.read_excel(r's',header=None).values
Dis0=np.load('s')

new_column=np.array([4 for i in range(len(p1))]) #安装高度4m
P=np.column_stack((p1,new_column))

#筛选北边的最外层
r_1=340
r_2=351
Ai=[]
for i in range(len(P)):
    if R_dis[i]>=r_1 and p1[i,1]>0:
        Ai.append(i)
        P[i,2]=P[i,2]+2
x_c=0
y_c=-250

Pi=[]
for i in range(len(P)):
    dis=sqrt((p1[i,0]-x_c)**2+(p1[i,1]-y_c)**2)
    if dis>=100:
        Pi.append(i)
PP=P[Pi]

np.save('sss',PP)
Dis0=np.load('ss')
D_0=[306]
ST_0=[9]

HR=3.5
H=4
L=W=6
H0=80
H1=8
delta_t=3 #步长

PW=F(L, W, H, D_0, ST_0, delta_t, H0, H1, HR, Dis0, PP, x_c, y_c)
print(np.mean(PW)*W*L*len(PP)/1000)

plt.scatter(p1[:,0],p1[:,1],c='b')
plt.scatter(p1[Ai,0],p1[Ai,1],c='r')

这段代码也是在计算WE问题的可用性Wp，但是有些不同。代码首先加载了附件，并定义了一些常量和变量，如板子的尺寸L、W、H，安装高度等。

然后，代码新增了一个安装高度的维度，将其添加到数组P中，用于标记每个点的安装高度，安装高度为4m。

接下来，代码筛选出北边最外层的点，并将这些点的安装高度加2，假设这些点高度更高，用于计算WE算法在北部的情况下的可用性。

然后，代码对于每个点进行判断，判断距离该点最远的板子与该点之间的距离是否小于100，如果小于100，则记录该点的信息，并筛选掉距离该点最远的板子与该点之间距离大于等于100的点。筛选后得到的数据存储为PP，并将其用于计算WE算法在北部最外层的情况下的可用性。

接下来，代码对PP和Dis0进行计算，并将计算结果存储为PW，最后输出计算结果。

最后，代码将所有点进行可视化，并将筛选出的最外层点标记为红色。

?

6、

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import random
from math import cos,sin,acos,asin,pi,exp,sqrt,floor
import time
from Function_3 import F

p1=pd.read_excel(r's',header=None).values
Dis0=np.load('s')

HR=3.5
H=4
L=W=6
H0=80
H1=8
delta_t=3 #步长

#导入每块板 在某时刻的 俯仰角与方位角
D_0=[306,337,0,31,34,92,122,153,184,214,245,275]
ST_0=[9,10.5,12,13.5,15]
x_c=-250
y_c=-150

m_n=floor(0.3*len(p1))
random_numbers=random.sample(list(range(len(p1))),m_n)
p11=p1[random_numbers]
for i in range(m_n):
    for j in range(m_n):
      Dis[i,j]=Dis0[random_numbers[i],random_numbers[j]]
start_time=time.time()
PW=F(L, W, H, D_0, ST_0, delta_t, H0, H1, HR, Dis, P11, x_c, y_c)

end_time=time.time()
elapsed_time_seconds=end_time-start_time
elapsed_time_minutes=round(elapsed_time_seconds/60,2)
print(f"代码运行时间:{elapsed_time_minutes}分钟")

print(PW)

这段代码首先通过pandas库读取了一个Excel文件，并将其转换为numpy数组p1。然后加载了之前计算出的Dis0距离矩阵，将其存储为数组Dis0。

代码中还定义了一些与WE算法相关的参数和变量，如板子的尺寸L、W、H，时间间隔delta_t，地球半径HR等。

接下来，代码对于p1中的数据进行筛选，将其中的部分数据随机取出来，生成新的数组p11，用于计算WE算法的可用性。

然后，代码调用了一个名为F的函数，用于计算WE算法在当前情况下的可用性。该函数对应于我们之前提到的WE算法的主要计算部分，通过输入一系列参数进行计算。

最后，代码输出了计算结果PW，并计算了代码运行时间。