计算机仿真试题和答案

控制系统计算机辅助设计综合实验指导实验名称：连续与离散系统校正实验，系统可控性与可观性实验，系统的simulink仿真实验陈茜编实验人：苏建聪学号：200830810122班级：08电气工程及其自动化1班信息工程系实验任务书1. 有一个单位负反馈控制系统，如果控制对象的传递函数为设计要求：① 相角裕度≥45°；② 当系统的输入信号是单位斜坡信号时，稳态误差ess ≤0.04。

③ 要求绘制出校正后系统和未校正系统的Bode 图及其闭环系统的单位阶跃响应曲线，并进行对比。

2. 有一个单位负反馈控制系统，如果控制对象的传递函数为：试设计一个串联滞后校正装置。

设计要求：①相角裕度≥45°；② 当系统的输入信号是单位斜坡信号时，稳态误差ess ≤0.04。

③ 要求绘制出校正后系统和未校正系统的Bode 图及其闭环系统的单位阶跃响应曲线，并进行对比。

3. 有一个单位负反馈控制系统，如果控制对象的传递函数为 ()()4+=s s k s G p试设计一个串联超前滞后校正装置，设计要求： ①相角裕度≥45°；② 当系统的输入信号是单位斜坡信号时，稳态误差ess ≤0.04。

③ 要求绘制出校正后系统和未校正系统的Bode 图及其闭环系统的单位阶跃响应曲线，并进行对比。

4. 系统结构图如图所示，其中，采样周期Ts=0.01s ，被控对象()()110+=s s s G ，()s G h 为零阶保持器。

用W 变换法设计一超前校正装置D(z)，使系统相位裕度γ≥50°，校验设计后系统的性能指标。

5. 系统结构图如图所示，其中，采样周期Ts=0.01s ，被控对象1)s(0.2s k )(+=s G ο，()s G h 为零阶保持器。

用对数频率法设计D(z)，使系统开环增益k ≥30(1/s)，截止频率ωc ≥15(1/s),相位裕度γ≥50 °1使 s 11se -1(s)-Tsh T G +≈=，求出未校正系统的开环系统的开环传递函数(s)(s)G G (s)0s =G ,的传递函数模型参数。

2作未校正系统的Bode 图，求出截止频率 和相位裕度 ，并与系统性能指标要求比较 3用串联超前校正，选择Wc=20(1/s)，确定D （s ）的参数。

4求出校正后系统的相位裕度γ，并检验是否满足要求，若不满足要求，则返回到3 5用双线性变换算法求出D （z ）. 6检查系统的性能6. 已知系统状态空间方程： u x x ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡+⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---=∙102011220001，[]x y 021=。

对系统进行可控性|、可观性分析以及极点配置控制器的设计。

实验要求：（1） 判别系统的可控性，求解系统的变换矩阵c Q 。

（2） 如果系统完全可控，导出系统的第一可控标准型。

（3） 判断系统的可观测性，求解系统的变换矩阵o Q （4） 如果系统完全可观测，导出系统的第一可观测标准型。

（5） 确定状态反馈矩阵K ，使系统的闭环极点配置在 []iip 3.24.33.24.32--+--=位置上，并绘制出配置后系统单位阶跃时间响应曲线。

7. 建立如图所示双闭环计算机控制直流调速系统的simulink 模型。

仿真观测转速和电流跟踪特性等性能指标，调整数字PI调节器参数，以满足系统的性能指标要求。

8. 有限拍随动系统的设计：1） 某一离散控制系统的被控队形的传递函数为：()()35+=s s s G p ，保持器为零阶保持器，采样周期,5.0s T =试设计单位速度输入时的有限拍控制器。

2） 某一离散控制系统的被控队形的传递函数为：()()35+=s s s G p ，保持器为零阶保持器，采样周期,5.0s T =试设计单位速度输入时的有限拍无纹波控制器。

9．用函数调用的方法，求出!!2!1n ⋅⋅⋅++，要求写出脚本式M 文件和函数式M 文件 10. 找出101——n 之间的满足（1）它是完全平方数； （2）在满足完全平方数的前提下，有两位数字相同。

统计满足条件的数的个数。

11 通过键盘输入一个整数，将该数变成按相反顺序的数，例如12345变成54321答案1.主程序:num=8000;den=conv([1,0],conv([1,4],[1,80]));G=tf(num,den); %未校正系统的开环传递函数[Gm,Pm,Wcg,Wcp]=margin(G); %1) 未校正系统的频域响应参数，计算需要相角欲度Pm(γ)w=0.1:0.1:10000; %确定频率的取值范围和频率采样的间隔值[mag,phase]=bode(G,w);magdb=20*log10(mag); %对数幅频、相频特性向量，计算需要magdb %(对数幅频响应值)phim1=45;deta=8; % 2）(4.1.9)式设置系统设计参数phim1( )；deta( )phim=phim1-Pm+deta; % (4.1.9)式相位超前角φbita=(1-sin(phim*pi/180))/(1+sin(phim*pi/180)); %求出β值n=find(magdb+10*log10(1/bita)<=0.0001); %找出满足magdb+10*log10(1/bita)<=0.0001%式的magdb向量所有的下标值wc=n(1); %通常magdb(1)+10*log10(1/bita)>0.0001;取的第1项为(比实际频率值大%10倍)，这是因为w=0.1:0.1:10000，而下标向量的值=1:1:100000w1=(wc/10)*sqrt(bita);w2=(wc/10)/sqrt(bita); % 3) (4.1.11)和(4.1.12)式numc=[1/w1,1];denc=[1/w2,1]; % (4.1.13)式，取K=1Gc=tf(numc,denc); %校正器的传递函数GmdB=20*log10(Gm); %4) 下面程序校验系统校正后的系统指标GcG=Gc*G;[Gmc,Pmc,wcgc,wcpc]=margin(GcG); %GcG是校正后系统的开环传递函数GmcdB=20*log10(Gmc);disp('未校正系统的开环传递函数和频域响应参数：h，γ，wc')G,[GmdB,Pm,Wcp],disp('校正装置传递函数和校正后系统开环传递函数')Gc,GcG,disp('校正后系统的频域响应参数：h，γ，wc')[GmcdB,Pmc,wcpc],disp('校正装置的参数T和β值：T，β')T=1/w1;[T,bita],bode(G,GcG);figure(2);margin(GcG);GcGc=feedback(GcG,1);step(numc,denc);运行结果：未校正系统的开环传递函数和频域响应参数：h，γ，wcTransfer function:8000--------------------s^3 + 84 s^2 + 320 sans =10.5268 15.8578 9.5715校正装置传递函数和校正后系统开环传递函数Transfer function:0.1458 s + 1-------------0.03602 s + 1Transfer function:1166 s + 8000-------------------------------------------0.03602 s^4 + 4.025 s^3 + 95.53 s^2 + 320 s校正后系统的频域响应参数：h，γ，wcans =16.0462 43.5208 13.7993校正装置的参数T和β值：T，βans =0.14580.2470截图：2.主程序:den=conv(den,[1,80]);G=tf(num,den); %未校正系统的传递函数margin(G);gamma_cas=45;delta=6;%设计要求的相角裕度，(4.1.21)式中的设置为gamma_l=gamma_cas+delta; %根据(4.1.21)式= +△w=0.01:0.01:1000; %设置计算频率点[mag,phase]=bode(G,w);n=find(180+phase-(gamma_l)<=0.1); %求出对应的频率，wgamma_l=n(1)/100; %n(1)的值与校正后实际频率相差100倍[mag,phase]=bode(G,wgamma_l); %求出在频率点的相位值rr=-20*log10(mag);beta=10^(rr/20); %根据(4.1.22)式：w2=wgamma_l/10;w1=beta*w2; %根据第④步，求出和numc=[1/w2,1];denc=[1/w1,1]; %根据(4.1.23)式，求出系统校正网络的传递函数Gc=tf(numc,denc)GcG=Gc*G %校正后系统的开环传递函数bode(G,GcG),figure(2),margin(GcG),betaGcGc=feedback(GcG,1);step(numc,denc);运行结果：Transfer function:3.344 s + 1-----------22.38 s + 1Transfer function:2.676e004 s + 8000---------------------------------------22.38 s^4 + 1881 s^3 + 7246 s^2 + 320 sbeta =0.1494截图：3.主程序:num=1600;den=conv([1,0],([1,4]));G=tf(num,den); %未校正系统的传递函数[h,gamma,wg,wc]=margin(G);h=20*log10(h); %未校正系统的频域响应参数w=0.001:0.001:100; %设置计算频率点[mag,phase]=bode(G,w); %未校正系统的幅频与相频特性值向量disp('未校正系统的参数：h,wc,γ');[ h,wc, gamma],gamma1=45;delta=6; %设计要求的相角裕度，(4.1.32)式中的设置为phim=gamma1-gamma+delta; % (4.1.32)式alpha=(1+sin(phim*pi/180))/(1-sin(phim*pi/180)); %(4.1.33)式magdb=20*log10(mag); %用分贝值表示的未校正系统的幅值裕度n=find(magdb+10*log10(alpha)<=0.0001); %找出满足magdb+10*log10(alpha)<=0.0001式%的magdb向量所有的下标值wc=n(1);wcc=wc/1000; % 的值与wc相差1000倍w3=wcc/sqrt(alpha);w4=sqrt(alpha)*wcc; %(4.1.34)式numc1=[1/w3,1];denc1=[1/w4,1]; %(4.1.35)式Gc1=tf(numc1,denc1); %超前校正部分的传递函数w1=wcc/10;w2=w1/alpha; %取= ；=numc2=[1/w1,1];denc2=[1/w2,1]; %(4.1.36)式Gc2=tf(numc2,denc2); %滞后校正部分的传递函数Gc12=Gc1*Gc2; %串联超前-滞后校正网络的传递函数GcG=Gc12*G; %校正后系统的传递函数[Gmc,Pmc,wcgc,wcpc]=margin(GcG);GmcdB=20*log10(Gmc);disp ('超前校正部分的传递函数'),Gc1,disp ('滞后校正部分的传递函数'),Gc2,disp ('串联超前-滞后校正网络的传递函数'),Gc12,disp('校正后系统的开环传递函数'),GcG,disp('校正后系统的性能参数：h,wc,γ及α值'),[GmcdB,wcpc,Pmc,alpha],bode(G,GcG),figure(2),margin(GcG),beta;GcGc=feedback(GcG,1);step(numc,denc);运行结果：未校正系统的参数：h,wc,γans =Inf 39.8991 5.7249超前校正部分的传递函数Transfer function:0.03902 s + 1--------------0.006604 s + 1滞后校正部分的传递函数Transfer function:0.1605 s + 1------------0.9484 s + 1串联超前-滞后校正网络的传递函数Transfer function:0.006263 s^2 + 0.1995 s + 1---------------------------0.006263 s^2 + 0.955 s + 1校正后系统的开环传递函数Transfer function:10.02 s^2 + 319.3 s + 1600----------------------------------------0.006263 s^4 + 0.98 s^3 + 4.82 s^2 + 4 s校正后系统的性能参数：h,wc,γ及α值ans =Inf 18.4720 25.6743 5.9083margin(sysw);[Gm,Pm,Wcg,Wcp]=margin(sysw);w=0.1:0.1:10000; %确定频率的取值范围和频率采样的间隔值[mag,phase]=bode(sysw,w);magdb=20*log10(mag); %对数幅频、相频特性向量，计算需要magdb%(对数幅频响应值)phim1=45;deta=8; % 2）(4.1.9)式设置系统设计参数phim1( )；deta( )phim=phim1-Pm+deta; % (4.1.9)式相位超前角φbita=(1-sin(phim*pi/180))/(1+sin(phim*pi/180)); %求出β值n=find(magdb+10*log10(1/bita)<=0.0001); %找出满足magdb+10*log10(1/bita)<=0.0001%式的magdb向量所有的下标值wc=n(1); %通常magdb(1)+10*log10(1/bita)>0.0001;取的第1项为(比实际频率值大%10倍)，这是因为w=0.1:0.1:10000，而下标向量的值=1:1:100000w1=(wc/10)*sqrt(bita);w2=(wc/10)/sqrt(bita); % 3) (4.1.11)和(4.1.12)式numc=[1/w1,1];denc=[1/w2,1]; % (4.1.13)式，取K=1disp('控制器W传递函数');Gc=tf(numc,denc); %校正器的传递函数disp('校正后系统的开环W传递函数');GcG=Gc*sysw %校正后系统的开环传递函数margin(GcG)sysc=c2d(Gc,0.01,'tustin')syszk=sysc*dsysg;margin(syszk)运行结果：Transfer function:4.983e-005 z + 4.967e-005-------------------------z^2 - 1.99 z + 0.99Sampling time: 0.01控制器W传递函数校正后系统的开环W传递函数Transfer function:-4.744e-008 s^3 - 0.005683 s^2 + 1.134 s + 1--------------------------------------------1.084 s^3 +2.084 s^2 + s + 1.232e-012Transfer function:1.05 z - 1.041--------------z - 0.9908Sampling time: 0.01截图：5. 主程序:s=0.01;[numh,denh]=pade(Ts/2,1);sysh=tf(numh,denh);numo=1;deno1=[0.2 1 0];sysgo1=tf(numo,deno1);sysgho=sysh*sysgo1;[kp,kv,ka,jienci]=wucha(sysgho);k=30/kv;sysgho=k*sysgho;margin(sysgho);[Gm,Pm,Wcg,Wcp]=margin(sysgho),h=20*log10(Gm);w=0.1:0.1:10000; %确定频率的取值范围和频率采样的间隔值[mag,phase]=bode(sysgho,w);magdb=20*log10(mag); %对数幅频、相频特性向量，计算需要magdb%(对数幅频响应值)phim1=50;deta=10; % 2）(4.1.9)式设置系统设计参数phim1( )；deta( )phim=phim1-Pm+deta; % (4.1.9)式相位超前角φbita=(1-sin(phim*pi/180))/(1+sin(phim*pi/180)); %求出β值n=find(magdb+10*log10(1/bita)<=0.0001); %找出满足magdb+10*log10(1/bita)<=0.0001%式的magdb向量所有的下标值wc=n(1); %通常magdb(1)+10*log10(1/bita)>0.0001;取的第1项为(比实际频率值大%10倍)，这是因为w=0.1:0.1:10000，而下标向量的值=1:1:100000w1=(wc/10)*sqrt(bita);w2=(wc/10)/sqrt(bita); % 3) (4.1.11)和(4.1.12)式numc=[1/w1,1];denc=[1/w2,1]; % (4.1.13)式，取K=1Gc=tf(numc,denc);disp('系统校正后的开环传递函数');GcG=Gc*sysghomargin(GcG);Gcd=c2d(Gc,0.01,'tustin')误差子程序:function [Kp,Kv,Ka,jienci]=wucha(sys) %sys：TF模型[num,den]=tfdata(sys,'v'); %求出模型的分子和分母的多项表达式[Kp,jienci]=jixian(num,den,0); %Kp：位置误差系数；jienci：系统类型[Kv,jienci]=jixian(num,den,1); %Kv：速度误差系数[Ka,jienci]=jixian(num,den,2); %Ka：加速度误差系数disp('系统误差系数')Kp,Kv,Ka,disp('系统类型是'),jienci极限子程序:function [lim,jienci]=jixian(knum,kden,nn)%求Kp、Kv、Ka，n=0，单位位置输入，n=1：单位速度输入，nn=2：单位加速度输入，lim：Kp、Kv、Ka，jienn：系统类型jien=length(kden);jienn=0; %jien：系统分母多项式的长度for i=1:jienif kden(jien-i+1)<=0.00001 %首先从多项式常数项开始判断，若为0，则系统类型值jienn=1+jienn; %jienn+1 jienn；然后依s升幂次序判断各项的系数elsebreak %若系数不为0，则退出循环，得到系统类型值endendjiemn=length(knum);knumk=knum(jiemn);jienci=jienn; %knumk：分子的常数项if nn==0&jienn==0 %单位位置输入，0型系统lim=knumk/kden(jien); %kden(jien)：分母的常数项else if nn==0lim='inf'; %单位位置输入，1、2型系统endendif nn==1&jienn==0 %单位速度输入，0型系统lim=0;else if nn==1&jienn==1 %单位速度输入，1型系统lim=knumk/kden(jien-1);else if nn==1 %单位速度输入，2型以上系统lim='inf';endendif nn==2&jienn<2 %单位加速度输入，0，1型系统lim=0;lim=0;else if nn==2&jienn==2 %单位加速度输入，2型系统lim=knumk/kden(jien-2);else if nn==2 %单位加速度输入，3型系统以上lim='inf';endendend运行结果：系统误差系数Kp =infKv =1Ka =系统类型是jienci =1Gm =6.7079Pm =19.6906Wcg =31.5241Wcp =11.7478系统校正后的开环传递函数Transfer function:-3.661 s^2 + 1434 s + 12000--------------------------------------------0.005232 s^4 + 2.319 s^3 + 91.46 s^2 + 400 sTransfer function:4.076 z - 3.755---------------z - 0.67916. 主程序:(1) A=[1 0 0 ;0 2 -2;-1 -1 0];B=[2;0;1];C=[1 2 0]; %系统状态方程模型n=length(A); %求系统阶次nQ=ctrb(A,B); %求解系统可变换矩阵m=rank(Q); %求系统可控性矩阵的秩mif m==n %if-else-end程序判断系统是否完全可控，满足m=n系统状态完全可控Ac1=inv(Q)*A*Q; %Bc1=inv(Q)*B;Cc1=C*Q; % ；disp('System is controllable.');disp('System First Controllable Canonnical Form is:'); Ac1,Bc1,Cc1disp('The Transformation Martrix is:');Qelse %m<n系统状态不完全可控disp('system state V ariable cannot be totally controlled');disp('The rank of System Controllable Martix is:');m %可控的状态变量数end(2) A=[1 0 0 ;0 2 -2;-1 -1 0];B=[2;0;1];C=[1 2 0]; %系统状态方程模型n=length(A); %求系统的秩nQ=obsv(A,C); %求解系统状态可观性矩阵,（6.2.3）式m=rank(Q); %求系统可观性矩阵的秩mif m==n %若m＝n则系统完全可观P=inv(Q); Ao1=inv(P)*A*P; Bo1=inv(P)*B;Co1=C*P; %（6.2.5）式disp('System is Observable.'); %系统是可观测的disp('System First Observable Canonnical Form is:'); Ao1, Bo1, Co1disp('The Transformation Martrix is:');Pelsedisp('system state V ariable cannot be totally Observed');disp('The rank of System Observable Martix is:'); mend(3) A=[1 0 0 ;0 2 -2;-1 -1 0];B=[2;0;1];C=[1 2 0]; %系统状态方程模型D=[0];disp('原系统的极点为');p=eig(A)'P=[-2;-3.4+sqrt(-5.29);-3.4-sqrt(-5.29)]; %P为所期望的极点K=place(A,B,P) %求状态反馈增益；或用acker(A,B,P)命令disp('配置后系统的极点为');p=eig(A-B*K)'disp('极点配置后的闭环系统为')sysnew=ss(A-B*K,B,C,D) %（6.3.1）式step(sysnew/dcgain(sysnew)) %绘制时间响应曲线运行结果：(1) System is controllable.System First Controllable Canonnical Form is:Ac1 =0 0 -21 0 00 1 3Bc1 =1Cc1 =2 -2 2The Transformation Martrix is:Q =2 2 20 -2 01 -2 0(2) System is Observable.System First Observable Canonnical Form is:Ao1 =0.0000 1.0000 0-0.0000 -0.0000 1.0000-2.0000 -0.0000 3.0000Bo1 =2.0000-2.00002.0000Co1 =1.0000 0 0The Transformation Martrix is:P =-2.0000 -2.0000 1.00001.5000 1.0000 -0.50001.0000 0.2500 -0.2500(3) 原系统的极点为p =-0.7321 2.7321 1.0000K =114.6250 299.1500 -217.4500配置后系统的极点为p =-3.4000 - 2.3000i -3.4000 + 2.3000i -2.0000 极点配置后的闭环系统为a =x1 x2 x3 x1 -228.3 -598.3 434.9x2 0 2 -2x3 -115.6 -300.2 217.5b =u1x1 2x2 0x3 1c =x1 x2 x3y1 1 2 0d =u1y1 0Continuous-time model.截图：7.截图：8. 主程序:(1)syms s T K z t kGs=(1/s)*K/s/(s+3); %ft=ilaplace(Gs);ftt=subs(ft,t,k*T); %求Gs的反变换，GoGpZ=(1-z^-1)*ztrans(ftt);GoGpZ=simplify(GoGpZ); %求，并化简Ez=(1-z^-1)^2; %Dcz=(1-Ez)/(Ez*GoGpZ);Dcz=simplify(Dcz); %(7.4.4)式，并化简Dcz=subs(Dcz,T,0.5);Dcz=subs(Dcz,K,5); % ，[dnum,dden]=numden(Dcz); %取得的分子和分母的多项式dnum=sym2poly(dnum);dden=sym2poly(dden); %转化为多项式系数向量disp('单位速度输入时的有限拍控制器Dc(z)=');[z p k]=tf2zp(dnum,dden);Dcz=zpk(z,p,k)disp('系统脉冲传递函数D(z)=C(z)/R(z)');Dz=1-Ez;Dz=simplify(Dz)(2) syms s T K z tGs=(1/s)*K/s/(s+1);ft=ilaplace(Gs);ftt=subs(ft,t,t*T);GoGpZ=(1-z^-1)*ztrans(ftt);GoGpZ=simplify(GoGpZ);GoGpZ=subs(GoGpZ,T,0.5);GoGpZ=subs(GoGpZ,K,5);[dnum,dden]=numden(GoGpZ);dnum=sym2poly(dnum);dden=sym2poly(dden);disp('广义对象的Z传递函数GoGp(z)=');[z p k]=tf2zp(dnum,dden);GoGpZ=zpk(z,p,k,0.5)syms a0 a1 b0 b1 uDz=u*(1+0.847*u)*(a0+a1*u); %(7.4.4)式Ez=(1-u)^2*(b0+b1*u); %(7.4.5)式ling=Dz-1+Ez;ling=collect(ling) %令：ling =- (847/1000)*( - (2.694/3.411409))(847/1000*a1+b1)*u^3+(-2*b1+b0+a1+847/1000*a0)*u^2+(b1-2*b0+a0)*u+b0-1eq1='b0-1=0';eq2='b1-2*b0+a0=0';eq3='-2*b1+b0+a1+847/1000*a0=0';eq4='847/1000*a1+b1=0'; [a0 a1 b0 b1]=solve(eq1,eq2,eq3,eq4,'a0,a1,b0,b1');syms uDz=u*(1+0.847*u)*(1.3311-0.7897*u);Ez=(1-u)^2*(1+0.6689*u);GoGpZ=0.53265*u*(1+0.8467*u)/(1-u)/(1-0.6065*u);Dcz=Dz/(Ez*GoGpZ);Dcz=simplify(Dcz);Dcz=collect(Dcz)运行结果：(1) 单位速度输入时的有限拍控制器Dc(z)=Zero/pole/gain:4.9784 (s-0.5) (s-0.2231)-------------------------(s-1) (s+0.6115)系统脉冲传递函数D(z)=C(z)/R(z)Dz =(2*z-1)/z^2(2) 广义对象的Z传递函数GoGp(z)=Zero/pole/gain:0.53265 (z+0.8467)------------------(z-1) (z-0.6065)Sampling time: 0.5ling =(847/1000*a1+b1)*u^3+(-2*b1+b0+a1+847/1000*a0)*u^2+(b1-2*b0+a0)*u+b0-1ling =0.6689ans =(847/1000*a1+b1)*u^3+(-2*b1+b0+a1+847/1000*a0)*u^2+(b1-2*b0+a0)*u+b0-1Dcz =-100/10653*(1000+847*u)*(-13311+7897*u)*(-2000+1213*u)/(-1+u)/(10000+6689*u)/(10000+8 467*u)9．主程序:n=input('n=');%输入一个n的值.s=0;%s赋初值0for i=1:ns=s+jiecheng(i);%调用jiecheng子程序ends阶乘子程序:function p=jiecheng(n);p=1;for i=1:np=p*i;endend运行结果：n=10s =4037913主程序:n=input('n=');%输入一个n的值.s=0;q=0;%s赋初值0p=1;for i=1:np=1;for q=1:ip=p*q;ends=s+p;ends运行结果：sn=10s =403791310.主程序:n=input('n=');%输入一个101至999之间的数字a=101;%对a赋初值101b=sqrt(n);%将n的开方值赋予bb=fix(b);%对b的值取整p=0;q=0;while a<n+1c=sqrt(a);%将a的开方值赋予cfor d=1:bif c==dp=p+1;t=a;x=rem(t,10);t=t/10;t=fix(t);y=rem(t,10);t=t/10;t=fix(t);if or(or(x==y,x==t),y==t)q=q+1;endendenda=a+1;endp%101-999之间完全平方数的个数q%101-999之前完全平方数且有两个数相同的个数运行结果：n=999p =21q =811主程序:x=input('x=');b=0;a=0;while x>0a=x;a=rem(a,10);x=x/10;x=fix(x);b=(b*10+a); endb运行结果：x=123456b =654321。