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镜坐标系?
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}分钟")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矩阵，即为最终的计算结果。
?
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，即为计算结果。
最后代码输出了程序运行的时间。
?
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。
最后，代码将生成的所有点进行可视化。
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，最后输出计算结果。
最后，代码将所有点进行可视化，并将筛选出的最外层点标记为红色。
?
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，并计算了代码运行时间。