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,并计算了代码运行时间。