import numpy as np from itertools import product class Lorry: def
__init__(self, station1, station2): # station1: P->D 换电站位置， station2：D->P 换电站位置
self.bet = 100 # 当前电池电量 self.p = 20 - station2 # 当前位置 self.t = 0 # 当前时间
self.station1 = station1 self.station2 = station2 self.n_charge = 0 # 充电次数 def
__str__(self): return "Lorry is at position {}, with battery {}, at time
{}".format(self.p, self.bet, self.t) def move(self, dt): #
0(P)___10(D)___20(P') # dt时间后的演化结果 if 0 <= self.p < 10 or np.isclose(self.p,
0): # P->D self.bet -= 1/2 * dt # 电池消耗 elif 10 <= self.p < 20 or
np.isclose(self.p, 10): # D->P self.bet -= 1/3 * dt # 电池消耗 if self.bet < 0 or
np.isclose(self.bet, 0): raise Exception("Lorry is out of battery") self.t =
self.t + dt # 时间变化 self.p = (self.p + dt * 1) % 20 # 位置变化 if np.isclose(self.p,
10) or np.isclose(self.p, 0): # 装卸货时间 self.t += 1 self.recharge() def
recharge(self): if (np.isclose(self.p, self.station1) or np.isclose(self.p, 20
- self.station2)) and 10 <= self.bet <= 12.6: # 逐渐降低阈值，测试是否存在可行解 self.bet = 100
# print("Lorry is recharged at time {}".format(self.t)) # gap.append(self.t) #
print(self) self.t += 2 self.n_charge += 1 # for (sta1, sta2) in
product(np.arange(0, 10, 1), np.arange(0, 10, 1)): for sta1 in np.arange(0.1,
10, 0.1): sta2 = sta1 lorry = Lorry(sta1, sta2) try: for T in
range(1000*60*10): lorry.move(dt=0.1) print('({:.1f}, {:.1f}) worked with
charge {:d} times'.format(sta1, sta2, lorry.n_charge)) except Exception as e:
print((sta1, sta2), e, end='\r') # # gap = [] # (sta1, sta2) = (5, 5) # lorry =
Lorry(sta1, sta2) # for T in range(1000*60*10): # lorry.move(dt=0.1) # #
print(lorry.n_charge) # # gap = [gap[i] - gap[i-1] for i in range(1, len(gap))]
# # gap = np.array(gap) # # print(gap.min()) # n_lorry = 75

load ('data.mat'); %读取数据 X=data; z=zscore(X) %数据标准化 M=cov(z) %协方差
[V,D]=eig(M); %求出协方差矩阵的特征向量、特征根 d=diag(D); %取出特征根矩阵列向量（提取出每一主成分的贡献率）
eig1=sort(d,'descend') %将贡献率按从大到小元素排列 v=fliplr(V) %依照D重新排列特征向量 S=0; i=0; while
S/sum(eig1)<0.85 i=i+1; S=S+eig1(i); end %求出累积贡献率大于85%的主成分 NEW=z*v(:,1:i)
%输出产生的新坐标下的数据 W=100*eig1/sum(eig1) figure(1) pareto(W); %画出贡献率的直方图
title('贡献率直方图'); %%%%%%因子分析%%%%%%%% %载荷 L =(v') * diag(sqrt(eig1)); %方差贡献 Var =
sum((L.^2)); %排序 Temp = [Var;1:length(Var)]; VarTemp =
flipud(sortrows(Temp',1)); L = L(:,VarTemp(:,2)); Var = VarTemp(:,1); %求维数的最小值
j = 0; Sum = 0; while Sum/sum(Var) <0.85 j= j +1; Sum = Sum + Var(j); end
L_main = L(:,1:j); %共同度 H = sum((L_main.^2),2); %方差贡献度 VarContr = 100 *
Var/sum(Var); figure(2) pareto(VarContr); title('因子旋转前方差贡献度');
[L_New,psi,T,stats,F]=factoran(z,j,'rotate','varimax','scores','Thomson');
%新的共同度 H = sum((L_New.^2),2); %新的方差贡献度 Var_New =sum(L_New.^2); VarContrNew =
100 * Var_New/sum(Var_New); figure(3) pareto(VarContrNew); title('因子旋转后方差贡献度');

50号与51号车间隔变大.

GitHub

Gitee