控制工程中的程序设计.docx
《控制工程中的程序设计.docx》由会员分享，可在线阅读，更多相关《控制工程中的程序设计.docx（13页珍藏版）》请在冰点文库上搜索。
控制工程中的程序设计
《控制工程中的程序设计》
大作业
班级：
学号：
姓名：
2015年5月1日
1.（20分）考虑以下迭代公式：
其中a、b为正数，要求：
1）编写程序求迭代的结果，迭代的终止条件为|xn+1-xn|≤10-5，迭代初值x0=1.0，迭代次数不超过500次。
2）如果迭代过程收敛于r，那么r的准确值是
，当（a,b）的值取（1,1）、（8,3）、（10,0.1）时，分别对迭代结果和准确值进行比较。
a=input（'a='）;%输入a的值
b=input（'b='）;%输入b的值
x=1.0;%迭代变量,确定初值
n=0;%初始迭代次数
whileabs（x-a/（b+x））>1e-5%终止条件
ifn<=500%迭代次数限制
x=a/（b+x）;
n=n+1;
end
end
n
x
r
（1）=（-b+sqrt（b*b+4*a））/2
r
（2）=（-b-sqrt（b*b+4*a））/2%输出r的准确值
s=r-x%迭代结果和准确值进行比较
实验结果：
1.
a=1
b=1
n=
12
x=
0.6180
r=
0.6180-1.6180
r=
0.6180-1.6180
s=
-0.0000-2.2361
2.
a=8
b=3
n=
12
x=
1.7016
r=
1.7016-1.6180
r=
1.7016-4.7016
s=
0.0000-6.4031
3.
a=10
b=0.1
n=
423
x=
3.1127
r=
3.1127-4.7016
r=
3.1127-3.2127
s=
-0.0000-6.3254
2.（15分）一系统可用下列方程组来表示：
从键盘输入m1、m2和θ的值，求a1、a2、N1和N2的值。
其中g取9.8，输入θ时以角度为单位。
要求：
定义一个求解线性方程组AX=B的函数文件，然后在命令文件中调用该函数文件。
函数fc.M文件:
functionX=fc（A,B）
%fcfc是求解线性方程的函数
%AA是未知矩阵的系数矩阵
X=A\B；
命令M文件：
clc;
m1=input（'m1='）;%输入m1的值
m2=input（'m2='）;%输入m2的值
theta=input（'theta='）;输入theta的值
x=theta*pi/180;
g=9.8
A=[m1*cos（x）-m1-sin（x）0
m1*sin（x）0cos（x）0
0m2-sin（x）0
00-cos（x）1];%输入矩阵A
B=[0;m1*g;0;m2*g];%输入矩阵B
X=fc（A,B）
实验结果：
m1=1
m2=1
theta=30
X=
7.8400
3.3948
6.7896
15.6800
3.（20分）绘制极坐标曲线ρ=a*sin（b+nθ），并分析参数a、b、n对曲线形状的影响。
theta=0:
pi/100:
2*pi;
a=input（'输入a='）;%输入a的值
b=input（'输入b='）;%输入b的值
n=input（'输入n='）;%输入n的值
rho=a*sin（b+n*theta）;%极半径
polar（theta,rho,'r'）%绘制二维极坐标图形。
采用控制变量法的办法，固定两个参数，变动第三个参数观察输出图象的变化。
结果：
a=b=n=2a=b=n=3
a=b=c=4
a=2,b=3,c=2a=2,b=4,c=2
a=2,b=5,c=2
a=3,b=2,c=2a=4,b=2,c=2
a=5,b=2,c=2
分析结果：
由这9个图知道，
当a,n固定时，图形的形状也就固定了，b只影响图形的旋转的角度;
当a,b固定时，n只影响图形的扇形数，特别地，当n是奇数时，扇叶数就是n,当是偶数时，扇叶数则是2n个;
当b,n固定时，a影响的是图形大小，特别地，当a是整数时，图形半径大小就是a。
4.（15分）若一个数等于它的各个真因子之和，则称该数为完数，如
6=1+2+3，所以6是完数。
求[1,500]之间的全部完数。
form=1:
500
sum=0;%sum初始化
forn=1:
m/2% 循环从1开始，步进是1，到m/2结束
ifrem（m,n）==0%无符号取余
sum=sum+n;
end
end
ifsum==m%判断如果sum和m相等
forn=1:
m/2%循环从1开始，步进是1，到m/2结束
ifrem（m,n）==0
sum=sum+n;
disp（['约数：
',num2str（n）]）;%显示约数
end
end
disp（['结果************',num2str（m）]）;%显示完数
end
end
运行结果：
约数：
1
约数：
2
约数：
3
结果************6
约数：
1
约数：
2
约数：
4
约数：
7
约数：
14
结果************28
约数：
1
约数：
2
约数：
4
约数：
8
约数：
16
约数：
31
约数：
62
约数：
124
约数：
248
结果************496
5.（15分）求非齐次线性方程组的通解。
A=[357;612;-8913]
B=[1-46;294;73-2]
C=[AB]%生成3*6矩阵
D=[A;B]%生成6*3矩阵
P=[C（2,1）C（2,4）C（2,5）;C（3,1）C（3,4）C（3,5）]%提取C的第2和3行,1/4/5列
AA=[2731;3522;9417]%系数矩阵
b=[6;4;2]
[m,n]=size（AA）
R=rank（AA）%求系数矩阵的秩
BB=[AAb]%增广矩阵
Rr=rank（BB）%求增广矩阵的秩
formatrat
ifR==Rr&R==n%n为未知数的个数，判断是否有唯一解
x=AA\b%唯一解
elseifR==Rr&Rx=AA\b%求特解
C=null（AA,'r'）%求AX=0的基础解系,所得C为n-R列矩阵，这n-R列即为对应的基础解系
%这种情形方程组通解xx=k（p）*C（:
P）（p=1n-R）
elseX='Nosolution!
'%判断是否无解
end
X='K1*C（:
1）+K2*C（:
2）+x'%非齐次线性方程组的通解
6.（15分）已知空调机组部分的数学模型为:
设采样时间T=1s，给定输入为单位阶跃信号，试采用PID控制方案。
y1（k）=0.848y1（k-1）+1.4687u1（k-1）+0.2381u2（k-1）
解：
clear all;
close all;
u1_1=0.0;u1_2=0.0;
u2_1=0.0;u2_2=0.0;
y1_1=0;y2_1=0;
x1=[0;0];x2=[0;0];x3=[0;0];
kp=0.285;
ki=0.098;
kd=0.0022;
error_1=[0;0];
ts=1;
for k=1:
1:
150
time（k）=k*ts;
R=[0;1];%PID Decouple Controller
u1（k）=kp*x1
（1）+kd*x2
（1）+ki*x3
（1）;
u2（k）=kp*x1
（2）+kd*x2
（2）+ki*x3
（2）;
u=[u1（k）,u2（k）];
if u1（k）>=10
u1（k）=10;
end
if u2（k）>=10
u2（k）=10;
end
if u1（k）<=-10
u1（k）=-10;
end
if u2（k）<=-10
u2（k）=-10;
end %Coupling Plant
yout1（k）= （0.848*y1_1+1.4687*u1_1+0.2381*u2_1）;
yout2（k）= （0.7678*y2_1+0.8536*u1_1+1.021*u2_1）;
error1（k）=R
（1）-yout1（k）;
error2（k）=R
（2）-yout2（k）;
error=[error1（k）;error2（k）];
u1_2=u1_1;u1_1=u
（1）;
u2_2=u2_1;u2_1=u
（2）;
y1_1=yout1（k）;
y2_1=yout2（k）;
x1=error; %Calculating P
x2=（error-error_1）/ts; %Calculating D
x3=x3+error*ts; %Calculating I
error_1=error;
end
figure
（1）;
plot（time,R
（1）,'k',time,yout1,'b'）;
hold on;
plot（time,R
（2）,'k',time,yout2,'r'）;
xlabel（'time（s）'）;ylabel（'rin,yout'）;
升级会员 