✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，
代码获取、论文复现及科研仿真合作可私信。
🍎个人主页：Matlab科研工作室
🍊个人信条：格物致知。
更多Matlab完整代码及仿真定制内容点击👇
🔥 内容介绍
随着电动汽车的普及，如何高效地进行充电成为了一个亟待解决的问题。在多充电站的情况下，如何实现电动汽车的有序充电，提高充电效率，降低充电成本，成为了一个重要的研究方向。本文将介绍基于粒子群算法实现IEEE33环境下的多充电站电动汽车群有序充电优化算法流程。
1.问题描述
考虑在IEEE33环境下，有多个充电站和多辆电动汽车。每个充电站都有一定数量的充电桩，每个充电桩的充电速度不同。每辆电动汽车都有一个充电需求，需要在规定时间内完成充电。同时，每个充电站的电力供应有限，需要在保证充电效率的前提下，合理分配电力资源。因此，需要设计一个有序充电方案，使得每辆电动汽车都能在规定时间内完成充电，同时最大化充电效率，最小化充电成本。
2.粒子群算法
粒子群算法是一种群体智能算法，模拟鸟群或鱼群等生物的集体行为，通过模拟群体的协作和竞争，寻找最优解。粒子群算法的基本思想是将每个解看作一个粒子，通过不断迭代，逐步接近最优解。粒子群算法具有收敛速度快、易于实现等特点，适用于解决复杂的优化问题。
3.算法流程
(1)初始化粒子群，包括每辆电动汽车的充电需求、每个充电站的充电桩数量和充电速度等信息。
(2)计算每个充电站的电力供应情况，并根据充电桩的充电速度和电动汽车的充电需求，确定每辆电动汽车的充电站和充电桩。
(3)根据每个充电站的电力供应情况，计算每辆电动汽车的充电时间和充电成本。
(4)根据充电时间和充电成本，更新每个粒子的适应度。
(5)根据粒子的适应度，更新粒子的速度和位置。
(6)重复步骤(2)-(5)，直到达到预设迭代次数或达到最优解。
4.实验结果
通过对IEEE33环境下的多充电站电动汽车群有序充电优化算法进行实验，得到了如下结果：
(1)在充电效率和充电成本方面，粒子群算法优于其他算法。
(2)随着粒子数的增加，算法的收敛速度更快，但计算时间也会增加。
(3)在充电桩数量和充电速度等参数的变化下，算法的性能表现稳定。
5.总结
本文介绍了基于粒子群算法实现IEEE33环境下的多充电站电动汽车群有序充电优化算法流程。该算法具有收敛速度快、易于实现等特点，在充电效率和充电成本方面优于其他算法。未来可以进一步优化算法，提高算法的性能表现。
📣 部分代码
%相关原始数据格式说明如下：%n——节点个数；n1——支路条数；isb——平衡节点号；H——PQ节点个数（为后面形成PVU存储PV节点初始电压用）；pr——误差精度。%B1——支路参数矩阵，其中第一列和第二列是起始节点编号和终点节点编号,第三列、第四列、第五列、第六列分别为：支路电阻、电抗、变压器变比、电纳。（不考虑电导）%B2——节点参数矩阵，其中第一列和第二列为节点编号和节点类型；第三列到第六列分别为：注入有功、注入无功、电压幅值、电压相位。%节点类型分类如下：“0”为平衡节点,“1”为PQ，“2”为PV节点；“3”为PQ（V）节点，“4”为PI节点。function [Ploss ,V_amp]=powerflow(car_1,car_2,car_3,P_flex)n=33 ;n=33 ;n1=32;isb=1;H=32; %%%%%%%%%%%%%18节点加DG PQV处理pr=0.01;v_amp=0;B1=[1 2 0.00922 0.0047i 1 0;2 3 0.00493 0.02511i 1 0;3 4 0.0366 0.01864i 1 0;4 5 0.03811 0.01941i 1 0;5 6 0.0819 0.0707i 1 0;6 7 0.01872 0.06188i 1 0;7 8 0.07114 0.02351i 1 0;8 9 0.103 0.074i 1 0;9 10 0.1044 0.074i 1 0;10 11 0.01966 0.0065i 1 0;11 12 0.03744 0.01238i 1 0;12 13 0.1468 0.1155i 1 0;13 14 0.05416 0.07129i 1 0;14 15 0.05910 0.0526i 1 0;15 16 0.07463 0.05450i 1 0;16 17 0.1289 0.1721i 1 0;17 18 0.0732 0.0574i 1 0;2 19 0.0164 0.01565i 1 0;19 20 0.15042 0.13554i 1 0;20 21 0.04095 0.04784i 1 0;21 22 0.07089 0.09373i 1 0;3 23 0.04512 0.03083i 1 0;23 24 0.08980 0.07091i 1 0;24 25 0.08960 0.07011i 1 0;6 26 0.0203 0.01034i 1 0;26 27 0.02842 0.01447i 1 0;27 28 0.1059 0.09337i 1 0;28 29 0.08042 0.07006i 1 0;29 30 0.05075 0.02585i 1 0;30 31 0.09744 0.0963i 1 0;31 32 0.03105 0.03619i 1 0;32 33 0.03410 0.05302i 1 0];B2=[1 0 0 0 1.05 0;2 1 -0.01 -0.006 1 0;3 1 -0.009 -0.004 1 0;4 1 -0.012 -0.008 1 0;5 1 -0.006 -0.003 1 0;6 1 -0.006 -0.002 1 0;7 1 -0.02 -0.01 1 0;8 1 -0.02 -0.01 1 0;9 1 -0.006 -0.002 1 0;10 1 -0.006 -0.0035 1 0;11 1 -0.0045 -0.003 1 0;12 1 -0.006 -0.0035 1 0;13 1 -0.006 -0.0035 1 0;14 1 -0.012 -0.008 1 0;15 1 -0.006 -0.001 1 0;16 1 -0.006 -0.002 1 0;17 1 -0.006 -0.002 1 0;18 1 -0.009 -0.004 1 0;19 1 -0.009 -0.004 1 0;20 1 -0.009 -0.004 1 0;21 1 -0.009 -0.004 1 0;22 1 -0.009 -0.004 1 0;23 1 -0.009 -0.005 1 0;24 1 -0.042 -0.02 1 0;25 1 -0.042 -0.02 1 0;26 1 -0.006 -0.0025 1 0;27 1 -0.006 -0.0025 1 0;28 1 -0.006 -0.002 1 0;29 1 -0.012 -0.007 1 0;30 1 -0.02 -0.06 1 0;31 1 -0.015 -0.007 1 0;32 1 -0.021 -0.01 1 0;33 1 -0.006 -0.004 1 0];for i=1:33B2(i,3)= B2(i,3)*P_flex; %负荷时变系数B2(i,4)= B2(i,4)*P_flex;endglobal positionB2(position(1),3)=B2(position(1),3)-car_1/10000; %电动汽车接入B2(position(1),4)=B2(position(1),4)-0.484*car_1/10000; %功率因素为0.9B2(position(2),3)=B2(position(2),3)-car_2/10000; %电动汽车接入B2(position(2),4)=B2(position(2),4)-0.484*car_2/10000; %功率因素为0.9B2(position(3),3)=B2(position(3),3)-car_3/10000; %电动汽车接入B2(position(3),4)=B2(position(3),4)-0.484*car_3/10000; %功率因素为0.9Y=zeros(n); %zeros就是生成一个全0的矩阵Times=1; %置迭代次数为初始值for i=1:n1p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/((B1(i,3)+B1(i,4))*B1(i,5));Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/(B1(i,3)+B1(i,4))+0.5*B1(i,6);%Y(q,q)=Y(q,q)+1/((B1(i,3)+B1(i,4))*B1(i,5)^2)+0.5*B1(i,6);%高压侧阻抗乘以变比平方 输入时注意低压侧在前end%disp('节点导纳矩阵：') ;Y;G=real(Y);B=imag(Y);OrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1); %将OrgS、DetaS初始化%创建OrgS，用于存储初始功率参数Q=0;PQV=0;x=1.655; %%%x1=6.7,x2=9.85,x=x1+x2；其中x1为定子漏抗，x2为转子漏抗xp=18.8; %%%xc=,xm=,xp=xc*xm/(xc-xm)；其中xc为机端并联电容器电抗，xm为激磁电抗h=0;for i=1:n %对PQ(V)节点的处理h=h+1;if i~=isb&&B2(i,2)==3Q(i)=-(B2(i,5))^2/xp+(-(B2(i,5))^2+sqrt((B2(i,5))^4-4*(B2(i,3))^2*x^2))/2*x;B2(i,4)=Q(i);B2(i,2)=1;PQV=h;endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Ig=0.01;Q=0;PI=0;h=0;for i=1:n %对PI节点的处理h=h+1;if i~=isb&&B2(i,2)==4Q(i)=sqrt(Ig^2*((B2(i,5))^2)-B(i,3)^2); %e=B2(i,5)，f=0，e^2+f^2=B2(i,5))^2，其中e和f为光伏发电系统接入节点电压的实部和虚部B2(i,4)=Q(i);B2(i,2)=1;PI=h;endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%h=0;j=0;for i=1:n %对PQ节点的处理if i~=isb&&B2(i,2)==1h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+B2(i,5)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))+B2(i,6)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5)); %Pi 书P57页11-45OrgS(2*h,1)=OrgS(2*h,1)+B2(i,6)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))-B2(i,5)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5)); %Qi 同上endendendfor i=1:n %对PV节点的处理if i~=isb&&B2(i,2)==2h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+ B2(i,5)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))+B2(i,6)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5)); %同上OrgS(2*h,1)=OrgS(2*h,1)+ B2(i,6)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))-B2(i,5)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5)); %同上endendendOrgS;%创建PVU 用于存储PV节点的初始电压PVU=zeros(n-H-1,2);t=0;for i=1:nif B2(i,2)==2t=t+1;PVU(t,1)=B2(i,5);PVU(t,2)=B2(i,6);endend%disp('PV节点初始值:电压、相位ei、fi：') ;PVU;%创建DetaS，用于存储有功功率、无功功率和电压幅值的不平衡量 书p58 11-46和11-47h=0;for i=1:n %对PQ节点的处理if i~=isb&&B2(i,2)==1h=h+1;DetaS(2*h-1,1)=B2(i,3)-OrgS(2*h-1,1);DetaS(2*h,1)=B2(i,4)-OrgS(2*h,1);endendt=0;for i=1:n %对PV节点的处理if i~=isb&&B2(i,2)==2h=h+1;t=t+1;DetaS(2*h-1,1)=B2(i,3)-OrgS(2*h-1,1);DetaS(2*h,1)=PVU(t,1)^2+PVU(t,2)^2-B2(i,5)^2-B2(i,6)^2;endend%disp('P、Q、V不平衡量：') ;DetaS;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%创建I，用于存储节点电流参数（计算雅克比矩阵用 具体见纸板公式推导）I=zeros(n-1,1);h=0;for i=1:nif i~=isbh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,5)+B2(i,6)*sqrt(-1));endendI;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%创建Jacbi(雅可比矩阵) 各元素是书中各元素求法乘个-1 所以△W=J△VJacbi=zeros(2*n-2);h=0;k=0;for i=1:n %对PQ节点的处理if B2(i,2)==1h=h+1;for j=1:nif j~=isbk=k+1;if i==j %对角元素的处理Jacbi(2*h-1,2*k-1)=-B(i,j)*B2(i,5)+G(i,j)*B2(i,6)+imag(I(h,1));Jacbi(2*h-1,2*k)=G(i,j)*B2(i,5)+B(i,j)*B2(i,6)+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));else %非对角元素的处理Jacbi(2*h-1,2*k-1)=-B(i,j)*B2(i,5)+G(i,j)*B2(i,6);Jacbi(2*h-1,2*k)=G(i,j)*B2(i,5)+B(i,j)*B2(i,6);Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k==(n-1) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行k=0;endendendendendk=0;for i=1:n %对PV节点的处理if B2(i,2)==2h=h+1;for j=1:nif j~=isbk=k+1;if i==j %对角元素的处理Jacbi(2*h-1,2*k-1)= -B(i,j)*B2(i,5)+G(i,j)*B2(i,6)+imag(I(h,1));Jacbi(2*h-1,2*k)= G(i,j)*B2(i,5)+B(i,j)*B2(i,6)+real(I(h,1));Jacbi(2*h,2*k-1)=2*B2(i,6);Jacbi(2*h,2*k)=2*B2(i,5);else %非对角元素的处理Jacbi(2*h-1,2*k-1)= -B(i,j)*B2(i,5)+G(i,j)*B2(i,6);Jacbi(2*h-1,2*k)= G(i,j)*B2(i,5)+B(i,j)*B2(i,6);Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k==(n-1) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行k=0;endendendendend%disp('雅克比矩阵为：') ;Jacbi; %列参数跟书上的相反 这里面第一列是对电压相位的偏导 第二列是对幅值的偏导%求解修正方程，获取节点电压的不平衡量DetaU=zeros(2*n-2,1);DetaU=Jacbi\DetaS;%disp('求解修正方程得：') ;DetaU;%修正节点电压j=0;for i=1:n %对PQ节点处理if B2(i,2)==1j=j+1;B2(i,5)=B2(i,5)+DetaU(2*j,1);B2(i,6)=B2(i,6)+DetaU(2*j-1,1);endendfor i=1:n %对PV节点的处理if B2(i,2)==2j=j+1;B2(i,5)=B2(i,5)+DetaU(2*j,1);B2(i,6)=B2(i,6)+DetaU(2*j-1,1);endendfor i=1:n %对PQ(V)节点的处理if i==PQVB2(i,2)=3;endendfor i=1:n %对PI节点的处理if i==PIB2(i,2)=4;endend%disp('修正后的B2阵：') ;B2;%开始循环**********************************************************************while max(abs(DetaU))>pr%先取绝对值在找出最大值来与误差精度比较for i=1:n %对PQ(V)节点的处理if i~=isb&&B2(i,2)==3Q(i)=-(B2(i,5))^2/xp+(-(B2(i,5))^2+sqrt((B2(i,5))^4-4*(B2(i,3))^2*x^2))/2*x;B2(i,4)=Q(i);endif B2(i,2)==3B2(i,2)=1;endendfor i=1:n %对PI节点的处理if i~=isb&&B2(i,2)==4Q(i)=sqrt(Ig^2*((B2(i,5))^2)-B(i,3)^2);B2(i,4)=Q(i);endif B2(i,2)==4B2(i,2)=1;endendOrgS=zeros(2*n-2,1);h=0;j=0;for i=1:n %对PQ节点的处理if i~=isb&&B2(i,2)==1h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+B2(i,5)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))+B2(i,6)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5)); %PiOrgS(2*h,1)=OrgS(2*h,1)+B2(i,6)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))-B2(i,5)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5)); %Qiendendendfor i=1:n %对PV节点的处理if i~=isb&&B2(i,2)==2h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+ B2(i,5)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))+B2(i,6)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5));OrgS(2*h,1)=OrgS(2*h,1)+ B2(i,6)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))-B2(i,5)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5));endendend%disp('修正后的迭代计算PQ、PV节点参数：') ;OrgS;%创建DetaSh=0;for i=1:n %对PQ节点的处理if i~=isb&&B2(i,2)==1h=h+1;DetaS(2*h-1,1)=B2(i,3)-OrgS(2*h-1,1);DetaS(2*h,1)=B2(i,4)-OrgS(2*h,1);endendt=0;for i=1:n %对PV节点的处理if i~=isb&&B2(i,2)==2h=h+1;t=t+1;DetaS(2*h-1,1)=B2(i,3)-OrgS(2*h-1,1);DetaS(2*h,1)=PVU(t,1)^2+PVU(t,2)^2-B2(i,5)^2-B2(i,6)^2;endend%disp('修正后的迭代计算PQ、PV节点不平衡量：') ;DetaS;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%创建II=zeros(n-1,1);h=0;for i=1:nif i~=isbh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,5)+B2(i,6)*sqrt(-1));endendI;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%创建JacbiJacbi=zeros(2*n-2);h=0;k=0;for i=1:n %对PQ节点的处理if B2(i,2)==1h=h+1;for j=1:nif j~=isbk=k+1;if i==j %对角元素的处理Jacbi(2*h-1,2*k-1)=-B(i,j)*B2(i,5)+G(i,j)*B2(i,6)+imag(I(h,1));Jacbi(2*h-1,2*k)=G(i,j)*B2(i,5)+B(i,j)*B2(i,6)+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));else %非对角元素的处理Jacbi(2*h-1,2*k-1)=-B(i,j)*B2(i,5)+G(i,j)*B2(i,6);Jacbi(2*h-1,2*k)=G(i,j)*B2(i,5)+B(i,j)*B2(i,6);Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k==(n-1) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行k=0;endendendendendk=0;for i=1:n %对PV节点的处理if B2(i,2)==2h=h+1;for j=1:nif j~=isbk=k+1;if i==j %对角元素的处理Jacbi(2*h-1,2*k-1)= -B(i,j)*B2(i,5)+G(i,j)*B2(i,6)+imag(I(h,1));Jacbi(2*h-1,2*k)= G(i,j)*B2(i,5)+B(i,j)*B2(i,6)+real(I(h,1));Jacbi(2*h,2*k-1)=2*B2(i,6);Jacbi(2*h,2*k)=2*B2(i,5);else %非对角元素的处理Jacbi(2*h-1,2*k-1)= -B(i,j)*B2(i,5)+G(i,j)*B2(i,6);Jacbi(2*h-1,2*k)= G(i,j)*B2(i,5)+B(i,j)*B2(i,6);Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k==(n-1) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行k=0;endendendendend%disp('修正后的雅克比矩阵：') ;Jacbi;DetaU=zeros(2*n-2,1);DetaU=Jacbi\DetaS;DetaU;%修正节点电压j=0;for i=1:n %对PQ节点处理if B2(i,2)==1j=j+1;B2(i,5)=B2(i,5)+DetaU(2*j,1);B2(i,6)=B2(i,6)+DetaU(2*j-1,1);endendfor i=1:n %对PV节点的处理if B2(i,2)==2j=j+1;B2(i,5)=B2(i,5)+DetaU(2*j,1);B2(i,6)=B2(i,6)+DetaU(2*j-1,1);endendfor i=1:n %对PQ(V)节点的处理if i==PQVB2(i,2)=3;endendfor i=1:n %对PI节点的处理if i==PIB2(i,2)=4;endendTimes=Times+1; %迭代次数加1enddisp('初始潮流计算结果如下：');disp('迭代次数：');Timesdisp('节点电压幅值如下（按节点由小到大的顺序）：');B2(:,5)%创建Sb，用于存储平衡节点功率Sb=zeros(1);for i=1:nif i==isbfor j=1:nSb=Sb+(B2(i,5)+sqrt(-1)*B2(i,6))*conj(Y(i,j))*conj(B2(j,5)+sqrt(-1)*B2(j,6));endendenddisp('初始平衡节点功率：')Sb%求解各条支路的功率Sij=zeros(n);for i=1:nfor j=1:nSij(i,j)=Sij(i,j)+(B2(i,5)+sqrt(-1)*B2(i,6))^2*conj(B1(1,4))+(B2(i,5)+sqrt(-1)*B2(i,6))*conj(Y(i,j))*(conj((B2(i,5)+sqrt(-1)*B2(i,6)))-conj((B2(j,5)+sqrt(-1)*B2(j,6))));endend%disp('线路功率分布：')Sij;Sij_P=real(Sij);Sij_Q=imag(Sij);for i=1:32P_flow(i)= abs( Sij_P(B1(i,1),B1(i,2)));endfor i=1:33for j=1:33if Sij_P(i,j)<0.0001Sij_P(i,j)=0;endif Sij_Q(i,j)<0.0001Sij_Q(i,j)=0;endendend%%%求解系统网损Ploss=0;for i=1:nfor j=1:nPloss=Ploss+B2(i,5)*B2(j,5)*(G(i,j)*cos(B2(i,6)-B2(j,6))+B(i,j)*sin(B2(i,6)-B2(j,6)));endendPloss=Ploss*10000;V_amp=B2(:,5);for i=1:n %电压罚函数V_33(i)=abs(1-V_amp(i)); %求与额定电压的差值endV_flu=sum(V_33);end
⛳️ 运行结果
🔗 参考文献
[1] 龚莉莉.基于利益链的配电网电动汽车充电优化策略研究[J].[2023-11-14].
[2] 李敏,苏小林,阎晓霞,等.多目标分层分区的电动汽车有序充放电优化控制[J].电网技术, 2015, 39(12):7.DOI:10.13335/j.1000-3673.pst.2015.12.033.
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
