资源描述
.\
附录： 学生实验报告表头格式
电力工程 学院学生实验报告
实验课程名称：
《电力系统计算机辅助分析》
开课实验室： 计算中心408 2014 年 12 月 11 日
年级、专业、班
学号
姓名
成绩
实验项目 名 称
电力系统计算机辅助分析
指 导 教 师
教 师 评 语
教师签名：
年 月 日
实验一 MATLAB 软件的基本操作及程序的基本结构
一、程序流程、源码及仿真结果
启动MATLAB软件，识别出常用的四个窗口――命令窗口（Command Window）、历史命令窗口（Command History）、工作空间浏览器（Workspace Browser）。熟悉Matlab工作环境。
2、上机求解下面的电路问题，脚本文件已经给出。
已知某个ＲＣ电路的端电压的表达式为：区间时，试绘制电压的波形。提示：在MATLAB的M文件编辑器中键入以下命令语句，并保存为ex.m：
％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％
Clear,clc,close
t=0:0.5:10;
u=6*exp(-2*t);
plot(t,u);title(‘RC电路的电压响应曲线’);xlabel(‘时间/s’);ylabel(‘电压/v’)
grid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5、 对于如下电路图，已知V，，，，，自己编写程序，画出和的波形（四个周期）。
提示：写出的向量表达式用于的计算（写出串联段和并联段的阻抗，利用分压公式），写出和的幅值和相角就能写出波形表达式，设定波形的长度（四个周波）和采样频率（两点间隔）后作图。注意画图时使用标题，横纵坐标及网格线等命令。
3.
%%%%%%%%%%%%%%%%%
clear,clc,close
mUsm=200; % 的幅值
aUsm=(pi/180)*(+90); % 的相角
Usm=mUsm*exp(j*aUsm); % 用于计算
w=100*pi;
r1=22;l=0.05;
r2=33;c=0.0001;
z1=r1+j*w*l;
z2=r2/(j*w*c*r2+1);
Ucm=Usm*z2/(z1+z2);
mUcm=abs(Ucm) %求幅值
aUcm=angle(Ucm) %求相角
T=4*2*pi/w; %用w表示 ，T表示四个周波
N=100; %画图的点数
dt=T/N; %两点间隔，即采样频率
t=0:dt:T;
for k=1:101 %画波形
us(k)=mUsm*cos(w*t(k)+aUsm);
uc(k)=mUcm*cos(w*t(k)+aUcm);
end
plot(t,us,t,uc) %注意区分两条曲线的线型和颜色
grid on
title()
xlabel(,fontsize,16) % fontsize’,16表示横坐标的字体大小为16
ylabel()
legend(us,uc)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
一、 讨论
1、 MATLAB软件是什么？有什么特点？
答： Matlab是一个高级的矩阵/阵列语言，它包含控制语句、函数、数据结构、输入和输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编写好一个较大的复杂的应用程（M文件）后再一起运行。
特点：1. 编程效率高，2. 用户使用方便，3. 扩充能力强，交互性好，4. 移植性和开放性很好，5. 语句简单，内涵丰富，6. 高效方便的矩阵和数组运算，7. 方便的绘图功能。
3. 如何获得命令或函数的帮助文档？
答：MATLAB 的各个函数，不管是内建函数、M文件函数、还是MEX文件函数等，一般它们都有M文件的使用帮助和函数功能说明，各个工具箱通常情况下也具有一个与工具箱名相间的M文件用来说明工具箱的构成内容等，在MATLAB命令窗口中，可以通过指令来获取这些纯文本的帮助信息。
通常能够起到帮助作用，获取帮助信息的指令有help、lookfor、which、doc、get、type等。
实验二 电力系统计算中常用的数值算法及电力网络的数学模型
1、 用因子表法求解电路方程中的节点电压；
A= [1/10+1/15 -1/10 -1/15 -1
-1/10 1/10+1/50+1/30 -1/50 0
-1/15 -1/50 1/15+1/50+1/30 0
1 0 0 0 ];
B=[0 0 0 10];
[n,m]=size(A);
for i=1:n
A(i,i)=1/A(i,i);
for j=i+1:n
A(i,j)=A(i,j)*A(i,i);
end
for k=i+1:n
for j=i+1:n
A(k,j)=A(k,j)-A(k,i)*A(i,j);
end
end
end
disp(矩阵A的因子表为：);
disp(A)
for i=1:n
B(i)=B(i)*A(i,i);
for j=i+1:n
B(j)=B(j)-A(j,i)*B(i);
end
end
for i=n-1:-1:1
for j=i+1:-1:2
B(j-1)=B(j-1)-A(j-1,i+1)*B(i+1);
end
end
disp(在因子表的基础上求解线性方程组的解为：x=);
disp(B)
2、将图2.1中的电压源分别改为2V、4V、6V、8V时，用因子表进行求解，体会该算法带来的益处；
① 电压源改为2V
②电压源改为4V
③电压源改为6V
④电压源改为8V
3、 上机验证课本例17.6－17.7的计算结果；
①例17.6
②例17.7
4、 设计系统变更时修改导纳矩阵的函数（参考课本P.323），并调试通过，给出导纳矩阵修正函数的调用规则及其主要的实现代码。
% n=input(请输入节点数:n=);
% nl=input(请输入支路数:nl=);
% B=input(请输入由支路参数形成的矩阵:B=);
% X=input(请输入由节点号及其对地阻抗形成的矩阵:X=);
n=4;nl=4;B=[1 2 0.08+0.4i 0 1 0;2 3 0.1+0.4i 0 1 0;
3 4 0.3i 0 1 0;1 3 0.12+0.5i 0 1 0];
X=[1 0;2 0;3 0;4 0]
Y=zeros(n);
for i=1:n
if X(i,2)~=0;
p=X(i,1);
Y(p,p)=1./X(i,2);
end
end
for i=1:nl
if B(i,6)==0
p=B(i,1); q=B(i,2);
else
p=B(i,2); q=B(i,1);
end
Y(p,q)=Y(p,q)-1./(B(i,3)*B(i,5));
Y(q,p)=Y(p,q);
Y(q,q)=Y(q,q)+1./(B(i,3)*B(i,5)^2)+B(i,4)./2;
Y(p,p)=Y(p,p)+1./B(i,3)+B(i,4)./2;
end
disp(导纳矩阵Y=:);
disp(Y)
B1=input(请输入追加支路参数形成的矩阵:B1=);
% B1=[3 4 0.3i 0 1.2 0];
% B1=[1 2 0.3i 0 1 0];
% B1=[1 4 0.3i 0 1 0];
% B1=[3 5 0.3i 0 1 0];
z=B1(3);
% if B1(6)==0
% p=B1(1); q=B1(2);
% else
% p=B1(2); q=B1(1);
% end
p=B1(1); q=B1(2);
if B1(6)==1
B1(5)=1/(B1(5));
end
if B1(2)>n
Y(p,p)=Y(p,p)+1/z; Y(p,q)=-1/z; Y(q,p)=-1/z; Y(q,q)=(1/z)+B1(4);
end
if Y(p,q)==0
Y(p,p)=Y(p,p)+1/z; Y(q,q)=Y(q,q)+1/z; Y(p,q)=-1/z; Y(q,p)=-1/z;
end
if Y(p,q)~=0&&B1(2)<=n
for i=1:nl
if B(i,1)==B1(1)&&B(i,2)==B1(2)
if B(i,5)==B1(5)
Y(p,p)=Y(p,p)+1/z; Y(q,q)=Y(q,q)+1/z;
Y(p,q)=Y(p,q)-1/z; Y(q,p)=Y(q,p)-1/z;
else
K1=B1(5);K=B(i,5);
Y(p,p)=Y(p,p)+(K1^2-K^2)*(1/z);Y(q,q)=Y(q,q); Y(p,q)=Y(p,q)-(K1-K)*(1/z); Y(q,p)=Y(q,p)-(K1-K)*(1/z);
end
end
end
end
if B1(2)<=n
n=n;
Y1=zeros(n);
for i=1:n
for j=1:n
Y1(i,j)=Y(i,j);
end
end
end
if B1(2)>n
n=n+1;
Y1=zeros(n);
for i=1:n
for j=1:n
Y1(i,j)=Y(i,j);
end
end
end
disp(新的导纳矩阵Y1=);
disp(Y1);
仿真结果：
X =
1 0
2 0
3 0
4 0
导纳矩阵Y=:
0.9346 - 4.2949i -0.4808 + 2.4038i -0.4539 + 1.8911i 0
-0.4808 + 2.4038i 1.0690 - 4.7568i -0.5882 + 2.3529i 0
-0.4539 + 1.8911i -0.5882 + 2.3529i 1.0421 - 7.5773i 0 + 3.3333i
0 0 0 + 3.3333i 0 - 3.3333i
请输入追加支路参数形成的矩阵:B1=[3 4 0.3i 0 1.2 0];
新的导纳矩阵Y1=
0.9346 - 4.2949i -0.4808 + 2.4038i -0.4539 + 1.8911i 0
-0.4808 + 2.4038i 1.0690 - 4.7568i -0.5882 + 2.3529i 0
-0.4539 + 1.8911i -0.5882 + 2.3529i 1.0421 - 9.0440i 0 + 4.0000i
0 0 0 + 4.0000i 0 - 3.3333i
思考题：
1、高斯消去法与因子表法有何异同？
答： 因子表法是以高斯消去法为基础求解线性方程组的一种方法，当线性方程组的系数矩阵不变，而常数项发生变化时，采用因子表法可避免系数矩阵的重复运算。因子表法首先要建立因子表，然后用因子表对常数项进行消去、回代，从而求出线性方程组的解。
2、节点导纳阵与节点阻抗阵之间有何关系，二者之间如何转换？
电力系统的网络方程式一般都用节点方程式表示，节点方程有导纳型和阻抗型两种。节点导纳矩阵具有对称、稀疏、可以根据网络接线图直接写出、易于修改的特点；节点阻抗矩阵是节点导纳矩阵的逆阵，它具有对称、满阵的特点，它不能由网络接线图直接写出，通常由支路追加法形成。
实验三 短路电流计算
合并课本第18章例程。即要求实现计算指定节点各种故障条件下，各支路的故障电流分布情况。其中对于实用算法，可只提供计算电抗。提示：尽量使用先前实验中已封装好的函数。
2.用新程序验证第18章的各例题。（注意短路电流曲线绘制时，至少包含5~7个点）
W=3;
switch W
case 1
clc
clear all
NF=1
for i1=1:NF
clear;
n1=3
n2=3
n0=3
nl1=5
nl2=5
nl0=5
Lf=1 %单相接地短路为1 横向故障为0
f=3
If=0
zf=0
zg=0
B1=[0 1 0.15i 0 1 0;0 2 0.075i 0 1 0;1 2 0.1i 0 1 0;1 3 0.1i 0 1 0;2 3 0.1i 0 1 0]
B2=[0 1 0.15i 0 1 0;0 2 0.075i 0 1 0;1 2 0.1i 0 1 0;1 3 0.1i 0 1 0;2 3 0.1i 0 1 0]
B0=[0 1 0.05i 0 1 0;0 2 0.025i 0 1 0;1 2 0.2i 0 1 0;1 3 0.2i 0 1 0;2 3 0.2i 0 1 0]
A1=[1 1 1;1 -1./2-i*sqrt(3)/2 -1./2+i*sqrt(3)/2;1 -1./2+i*sqrt(3)/2 -1./2-i*sqrt(3)/2]
Vcs=[1;1;1]
Z1=zeros(n1);Z2=zeros(n2);Z0=zeros(n0);Y1=zeros(n1);Y2=zeros(n2);Y0=zeros(n0);V1=zeros(n1,1);V2=zeros(n2,1);
V0=zeros(n0,1);I1=zeros(nl1,1);I2=zeros(nl2,1);I0=zeros(nl0,1);
for m1=1:3
m=0;
if m1==1
nl=nl1;B=B1;Z=Z1;n=n1;
elseif m1==2
nl=nl2;B=B2;Z=Z2;n=n2;
elseif m1==3
nl=nl0;B=B0;Z=Z0;n=n0;
end
m=0;Z=zeros(n);V=zeros(n);I=zeros(nl);
for k1=1:nl
p=B(k1,1);q=B(k1,2);
if B(k1,6)==0
k=1./B(k1,5);
else k=B(k1,5);
end
if p==0
if q>m %追加接地树支
Z(q,q)=B(k1,3);m=m+1;
else %追加接地连支
for i=1:m
Z(i,m+1)=-Z(i,q);Z(m+1,i)=-Z(q,i);
end
Z(m+1,m+1)=Z(q,q)+B(k1,3);
for i=1:m
for j=1:m
Z(i,j)=Z(i,j)-Z(i,m+1)*Z(m+1,j)./Z(m+1,m+1);
end
Z(i,m+1)=0
end
for i=1:m+1
Z(m+1,i)=0
end
end
else if q>m %追加不接地树支
for i=1:m
Z(i,q)=Z(i,p)*k;Z(q,i)=Z(p,i)*k;
end
Z(q,q)=k^2*Z(p,p)+k^2*B(k1,3);m=m+1;
else
for i=1:m %追加不接地连支
Z(i,m+1)=k*Z(i,p)-Z(i,q);Z(m+1,i)=k*Z(p,i)-Z(q,i);
end
Z(m+1,m+1)=k^2*Z(p,p)+Z(q,q)-2*k*Z(p,q)+k^2*B(k1,3);
for i=1:m
for j=1:m
Z(i,j)=Z(i,j)-Z(i,m+1)*Z(m+1,j)./Z(m+1,m+1);
end
Z(i,m+1)=0;
end
for i=1:m+1
Z(m+1,i)=0;
end
end
end
end
if m1==1
Z1=Z;
elseif m1==2
Z2=Z;
elseif m1==3
Z0=Z;
end
end %以上程序求出正,负,零序对应的阻抗阵
if If==0 %求出横向故障的阻抗阵
Z1(f,f)=Z1(f,f);Z2(f,f)=Z2(f,f);Z0(f,f)=Z0(f,f);
else %求出纵向故障的阻抗阵
Z1(f,f)=Z1(f,f)+Z1(1,1)-2*Z1(f,1);Z2(f,f)=Z2(f,f)+Z2(1,1)-2*Z2(f,1);Z0(f,f)=Z0(f,f)+Z0(1,1)-2*Z0(f,1);
end %-----------------故障点自阻抗--------------------------------
if Lf==1 %求出正负零序短路节点的电流标幺值
Z4=Z2(f,f)+Z0(f,f)+3*zf;k2=1;;k0=1;I1(f)=Vcs(f)./(Z1(f,f)+Z4);I2(f)=k2*I1(f);I0(f)=k0*I1(f);a1=I1(f);a2=I2(f);a0=I0(f);
elseif Lf==2 %两相短路
Z4=Z2(f,f)+2*zf;k2=-1;k0=0;
I1(f)=Vcs(f)./(Z1(f,f)+Z4);I2(f)=k2*I1(f);
I0(f)=k0*I1(f);a1=I1(f);
a2=I2(f);a0=I0(f);
elseif Lf==3 %两相接地短路
Z4=zf+(Z2(f,f)+zf)*(Z0(f,f)+zf+3*zg)./(Z2(f,f)+Z0(f,f)+2*zf+3*zg);k2=-(Z0(f,f)+zf+3*zg)./(Z2(f,f)+Z0(f,f)+2*zf+3*zg);
k0=-(Z2(f,f)+zf)./(Z2(f,f)+Z0(f,f)+2*zf+3*zg);I1(f)=Vcs(f)./(Z1(f,f)+Z4);I2(f)=k2*I1(f);I0(f)=k0*I1(f);a1=I1(f);
a2=I2(f);a0=I0(f);
elseif Lf==4 %单相断线
Z4=zf+(Z2(f,f)+zf)*(Z0(f,f)+zf)./(Z2(f,f)+Z0(f,f)+2*zf);k2=-(Z0(f,f)+zf)./(Z2(f,f)+Z0(f,f)+2*zf);
k0=-(Z2(f,f)+zf)./(Z2(f,f)+Z0(f,f)+2*zf);I1(f)=Vcs(f)./(Z1(f,f)+Z4);I2(f)=k2*I1(f);I0(f)=k0*I1(f);a1=I1(f);a2=I2(f);a0=I0(f);
elseif Lf==5 %两相断线
Z4=Z2(f,f)+Z0(f,f)+3*zf;k2=1;k0=1;I1(f)=Vcs(f)./(Z1(f,f)+Z4);I2(f)=k2*I1(f);I0(f)=k0*I1(f);a1=I1(f);a2=I2(f);a0=I0(f);
end
disp(故障点处的正序电流标幺值为:);
disp(a1);
disp(故障点处的负序电流标幺值为:);
disp(a2);
disp(故障点处的零序电流标幺值为:);
disp(a0);
for i=1:n1 %求出各节点的电压标幺值
V1(i)=Vcs(i)-I1(f)*Z1(i,f); %正序
end
for i=1:n2
V2(i)=-I2(f)*Z2(i,f); %负序
end
for i=1:n0
V0(i)=-I0(f)*Z0(i,f); %零序
end
disp(各节点的正序电压标幺值(节点号从小到大排):);
disp(V1);
disp(各节点的负序电压标幺值(节点号从小到大排):);
disp(V2);
disp(各节点的零序电压标幺值(节点号从小到大排):);
disp(V0);
for i=1:3 %求各支路的各序电流标幺值
if i==1
B=B1;V=V1;nl=nl1;
elseif i==2
B=B2;V=V2;nl=nl2;
elseif i==3
B=B0;V=V0;nl=nl0;
end
for j=1:nl
if B(j,6)==0
k=B(j,5);
else k=1./B(j,5);
end
p=B(j,1);q=B(j,2);
if p==0
e=0;b=B(j,3);
I(j)=(e-V(q)./k)./b; %接地支路电流
else
I(j)=(V(p)-V(q)./k)./B(j,3); %不接地支路电流
end
end
if i==1
I1=I;
elseif i==2
I2=I;
elseif i==3
I0=I;
end
end
disp(各支路的正序电流标幺值为(顺序同您输入B时一样):);
for i=1:nl
disp(I1(i));
end
disp(各支路的负序电流标幺值为(顺序同您输入B时一样):);
for i=1:nl
disp(I2(i));
end
disp(各支路的零序电流标幺值为(顺序同您输入B时一样):);
for i=1:nl
disp(I0(i));
end
for i=1:3 %求故障点A,B,C三相的总电流,各个节点A,B,C三相各相的总电压,各支路A,B,C三相各相的总电流
if i==1 %故障点电流
s0=a0;s1=a1;s2=a2;
elseif i==2 %各序支路电压
s0=V0;s1=V1;s2=V2;
else %各序支路电流
s0=I0;s1=I1;s2=I2;
end
for j=1:3 %不对称系数
CA=s0*A1(j,1)+s1*A1(j,2)+s2*A1(j,3);
if j==1
A=CA;
elseif j==2
B=CA;
elseif j==3
C=CA;
end
end
if i==1
Iadl=A;Ibdl=B;Icdl=C;
elseif i==2
Vajy=A;Vbjy=B;Vcjy=C;
else
Iazl=A;Ibzl=B;Iczl=C;
end
end
disp(故障点处A相电流Iadl为:);disp(Iadl);disp(故障点处B相电流Ibdl为:);disp(Ibdl);
disp(故障点处C相电流Icdl为:);disp(Icdl);disp(各节点处A相电压Vajy分别为(节点号从小到大排):);disp(Vajy);
disp(各节点处B相电压Vbjy分别为(节点号从小到大排):)disp(Vbjy);
disp(各节点处C相电压Vcjy分别为(节点号从小到大排):);disp(Vcjy);
disp(各支路A相电流Iazl分别为(顺序同你输入B时一样):);
for i=1:nl
disp(Iazl(i));
end
disp(各支路B相电流Ibzl分别为(顺序同你输入B时一样):);
for i=1:nl
disp(Ibzl(i));
end
disp(各支路C相电流Iczl分别为(顺序同你输入B时一样):);
for i=1:nl
disp(Iczl(i));
end
end
case 2
NF=1;
n=4;
nl=5;
B=[0 1 0.2i 0 1 0;0 2 4i 0 1 0;1 3 0.51i 0 1 0;2 3 0.59i 0 1 0;3 4 1.43i 0 1 0] ;D=[4,0];V0=[1;1;1;1];
m=0;Z=zeros(n);V=zeros(n);I=zeros(nl);
for k1=1:nl
p=B(k1,1);q=B(k1,2);
if B(k1,6)==0
k=1./B(k1,5);
else k=B(k1,5);
end
if p==0
if q>m
Z(q,q)=B(k1,3);m=m+1;
else
for i=1:m,
Z(i,m+1)=-Z(i,q);Z(m+1,i)=-Z(q,i);
end
Z(m+1,m+1)=Z(q,q)+B(k1,3);
for i=1:m
for j=1:m
Z(i,j)=Z(i,j)-Z(i,m+1)*Z(m+1,j)./Z(m+1,m+1);
end
Z(i,m+1)=0;
end
for i=1:m+1
Z(m+1,i)=0;
end
end
else if q>m
for i=1:m
Z(i,q)=Z(i,p)*k;Z(q,i)=Z(p,i)*k;
end
Z(q,q)=k^2*Z(p,p)+k^2*B(k1,3);
m=m+1;
else
for i=1:m
Z(i,m+1)=k*Z(i,p)-Z(i,q);Z(m+1,i)=k*Z(p,i)-Z(q,i);
end
Z(m+1,m+1)=k^2*Z(p,p)+Z(q,q)-2*k*Z(p,q)+k^2*B(k1,3);
for i=1:m
for j=1:m
Z(i,j)=Z(i,j)-Z(i,m+1)*Z(m+1,j)./Z(m+1,m+1);
end
Z(i,m+1)=0;
end
for i=1:m+1
Z(m+1,i)=0;
end
end
end
end
for k=1:NF
I(D(k,1),D(k,1))=V0(D(k,1),1)./(Z(D(k,1),D(k,1))+D(k,2));ft=num2str(D(k,1));ts1=(点短路时);ts2=(电流的标么值If=);
dn=strcat(ft,ts1,ts2);disp(dn);disp(I(D(k,1),D(k,1)));
for i=1:n
V(i,i)=V0(i,1)-I(D(k,1),D(k,1))*Z(i,D(k,1));
end
for i=1:nl
if B(i,6)==0
K=B(i,5);
else k=1./B(i,5);
end
p=B(i,1);q=B(i,2);
if p==0
e=0;b=B(i,3);
I(i,i)=(e-V(q,q)./k)./b;
else
I(i,i)=(e-V(p,p)-V(q,q)./k)./B(i,3);
end
end
disp(各节点的电压标么值U为：);
for i=1:n
disp(V(i,i));
end
disp(各支路短路电流的标么值I为：);
for i=1:nl
disp(I(i,i));
end
end
case 3
clc
NF=1;
n=3;
nl=5;
Sb=100;
D=[3 0];
vb=[115 115 115];
V0=[1;1;1];
B=[0 1 0.15i 0 1 0 100 1 1;0 2 0.075i 0 1 0 200 1 1;1 2 0.1i 0 1 0 0 0 2;1 3 0.1i 0 1 0 0 0 2;2 3 0.1i 0 1 0 0 0 2];
N=5;
m=0;Z=zeros(n);sq3=sqrt(3);s=zeros(N,1);
for k1=1:nl
p=B(k1,1);q=B(k1,2);
if B(k1,6)==0
k=1./B(k1,5);
else k=B(k1,5);
end
if p==0
if q>m
Z(q,q)=B(k1,3);m=m+1;
else
for i1=1:m,
Z(i1,m+1)=-Z(i1,q);Z(m+1,i1)=-Z(q,i1);
end
Z(m+1,m+1)=Z(q,q)+B(k1,3);
for i1=1:m
for j=1:m
Z(i1,j)=Z(i1,j)-Z(i1,m+1)*Z(m+1,j)./Z(m+1,m+1);
end
Z(i1,m+1)=0;
end
for i1=1:m+1
Z(m+1,i1)=0;
end
end
else if q>m
for i1=1:m
Z(i1,q)=Z(i1,p)*k;Z(q,i1)=Z(p,i1)*k;
end
Z(q,q)=k^2*Z(p,p)+k^2*B(k1,3);
m=m+1;
else
for i1=1:m
Z(i1,m+1)=k*Z(i1,p)-Z(i1,q);Z(m+1,i1)=k*Z(p,i1)-Z(q,i1);
end
Z(m+1,m+1)=k^2*Z(p,p)+Z(q,q)-2*k*Z(p,
展开阅读全文
相关搜索