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.5def 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*4yita_ref=0.92 #镜面反射率S=L*Wnew_column=np.array([H for i in range(len(p1))]) #安装高度4mP=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)**2b=0.5055+0.00595-(6.5-3)**2c=0.2711+0.01858*(2.5-3)**2DNI=1.366*(a+b*exp(-c/sin_as)) #得到法向直接辐射辐照度 DNIA0=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-Dinl=-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=0empty=0barr_tower=0barr_s=0barr_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.02for 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、zx0,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+=1continue#########二、判断入射光线是否被其他光镜遮挡##############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+=1if_barr=1break##################三、对应的反射###############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+=1if_barr=1break###########四、是否吸收################if if_barr==0:a,b,c=g_rx0,y0,z0=Di_ddelta_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+=1else:empty+=1#这里是这个时间点,计算第i个面板的数值,并记录,每个列表是1745的长度yita_sb=1-(barr_r+barr_s+barr_tower)/Nyita_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=1else:yita_trunc=(light)/(N-barr_s-barr_r-barr_tower)yita=yita_sb*yita_cos*yita_at*yita_trunc*yita_refYita_sb[sti,i]=yita_sbYita_cos[sti, i] = yita_cosYita_at[sti, i] = yita_atYita_trunc[sti, i] = yita_truncYita[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).Tp1=p1Dis0=np.load('D')xy_len=range(-250,251,50)delta_t=1EP=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:continueEP0=0Pi=[]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)/1000print(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:nW(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:nW(k,1)=W(k-1,1)+r;W(k,2)=W(k-1,2);k=k+1;
endn1=fix(2*R/(r*sqrt(3)/2))+1;
%从第三排开始
for i=3:n1if 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:nW(k,1)=W(k-1,1)+r;W(k,2)=W(k-1,2);k=k+1;endelse %偶数行%先写该行的第一个点W(k,1)=W2(1,1);W(k,2)=W2(1,2)-sqrt(3)*r*(i-2)/2;k=k+1;for j=1:nW(k,1)=W(k-1,1)+r;W(k,2)=W(k-1,2);k=k+1;endend
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^2a=a+1;W_e(a,:)=W(i,:);end
endscatter(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 Fp1=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=-250Pi=[]
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 Fp1=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=-150m_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,并计算了代码运行时间。