《经济控制论》课程研究报告姓名学号完成日期成绩联系方式一、引言：随着经济控制研究的深入，经济系统的建模、性能分析和优化设计等问题越来越趋向于纵深化和复杂化，MATLAB是具有很强数值计算、符号运算方针和图形显示功能的计算分析软件。

本次课程研究的目的，旨在加强对MATLAB的使用，学习使用MATLAB来解决经济控制系统的稳定性分析，状态观测器的设计以及经济控制系统的优化方法。

二、数据资料：1、宏观经济各个变量总量平衡统计数据表数据来源：根据中国统计年鉴整理。

更多数据请登录国家统计局网站查找: / 2、根据表中统计数据得到如下回归方程：消费函数 C(k) = 278.85 +0.5654Y(k-1) 投资函数 I(k) = 1338+3.3664[C(k)-C(k-1)]年度 国民生产总值（亿元） Y （k ）全社会固定资产投资(亿元) I(k)居民消费支出亿元）C(k)财政支出 （亿元）G （k ） 总量平衡 (I+C+G)/Y1986 10201.4 3120.6 5175.0 2204.91 1.0293 1987 11954.9 3791.7 5961.2 2362.18 1.0134 1988 14922.3 4753.8 7633.1 2491.21 0.9970 1989 16917.8 4410.4 8523.5 2823.78 0.9314 1990 18598.4 4517.0 9113.2 3083.59 0.8987 1991 21662.5 5594.5 10315.9 3386.62 0.8908 1992 26651.9 8080.1 12459.3 3742.20 0.9111 1993 34560.5 13072.3 15682.4 4642.30 0.9663 1994 46670.0 17042.1 20809.8 5792.62 0.9352 1995 57494.9 20019.3 26944.5 6823.72 0.9355 1996 66850.522913.532152.37937.55 0.9425政府购买 G(k) = 68.617k 2 - 271.51k + 2590，k 为年份，设1986年为0，依次递增。

三、以政府购买为控制变量，建立动态乘数－加速数经济模型的状态方程。

解答：系统的状态方程：输出方程：其中，a=3.3664;b=0.5654[C0 I0]=[278.85; 2.2767e+003]四、判定该经济系统的可控性、可观性、稳定性。

syms a bA=[b b;a*b-a a*b] B=[b ;a*b] C=[1 1] D=1%判状态能控性 P2=[B A*B] det(P2) rank(P2)[])()()(11)()()()()1()1(00k U k I k C k y I C k U ab b k I k C ab aab b b k I k C +⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡-=⎥⎦⎤⎢⎣⎡++if rank(P2)==2 disp('-------系统能控--------')else disp('--------系统不能控--------')end%判系统输出的能控性Q2=[C*B C*A*B D]rank(Q2)if rank(Q2)==1 disp('---------系统输出能控----------') else disp('---------系统输出不能控---------')end%判定系统的能观性R2=[C;C*A]rank(R2)if rank(R2)==2 disp('---------系统能观----------') else disp('---------系统不能观----------')end运行结果：A =[ b, b][ a*b - a, a*b]B = ba*bC =[ 1 1]D = 1P2 =[ b, a*b^2 + b^2][ a*b, a^2*b^2 - b*(a - a*b)]ans =-a*b^2ans =2-------系统能控--------Q2 =[ b + a*b, b*(b - a + a*b) + a*b*(b + a*b), 1] ans =1---------系统输出能控----------R2 =[ 1, 1][ b - a + a*b, b + a*b]ans =2---------系统能观----------%判定原系统的稳定性x=sym('x')L=sym('lambda')system_root=eig(A)pretty(system_root)if abs('system_root')<1 disp('------开环系统稳定------')else disp('------开环系统不稳定------')end运行结果：x =xL =lambdasystem_root =b/2 + (a*b)/2 - (b*(b - 4*a + 2*a*b + a^2*b))^(1/2)/2b/2 + (a*b)/2 + (b*(b - 4*a + 2*a*b + a^2*b))^(1/2)/2+- -+| 2 1/2 || b a b (b (b - 4 a + 2 a b + a b)) || - + --- - ------------------------------- || 2 2 2 || || 2 1/2 || b a b (b (b - 4 a + 2 a b + a b)) || - + --- + ------------------------------- || 2 2 2 |+- -+ ------开环系统不稳定------五、应用状态反馈改善系统性能。

%引入状态反馈syms f0 f1F=[f0 f1]L*eye(size(A)) -( A-B*F)poly(A-B*F)collect(poly(A-B*F),x)%判定稳定的区域a0=sym('a0')a1=sym('a1')a0=a*b-a*b*f1a1=b*f0-b-a*b+a*b*f1a2=1-a0a3=1+a1+a0a4=1-a1+a0A1= A-B*FP22=[B A1*B]rank(P22)if rank(P22)==2 disp('------闭环系统能控，能控性不变------') else disp('------闭环系统系统不能控------')enda=3.3664;b=0.5654;A=[b b;a*b-a a*b]B=[b ;a*b]C= [278.85; 2.2767e+003]D=1roots(poly(A))abs(roots(poly(A)))if abs(roots(poly(A)))<1disp('------开环系统稳定------')else disp('------开环系统不稳定------')endP2=[B A*B]rank(P2)if rank(P2)==2 disp('------开环系统能控------')else disp('------开环系统不能控------')end%系统稳定区域的判定a0=a*b-a*b*f1a1=b*f0-b-a*b+a*b*f1'闭环系统稳定条件: 1-a0>0,1+a1+a0>0 ,1-a1+a0>0'a2=1-a0a3=1+a1+a0a4=1-a1+a0运行结果：F =[ f0, f1]ans =[ lambda - b + b*f0, b*f1 - b][ a - a*b + a*b*f0, lambda - a*b + a*b*f1]ans =x^2 + (b*f0 - a*b - b + a*b*f1)*x + (a*b - a*b*f1)*(b - b*f0) + (b - b*f1)*(a - a*b + a*b*f0)ans =x^2 + (b*f0 - a*b - b + a*b*f1)*x + (a*b - a*b*f1)*(b - b*f0) + (b - b*f1)*(a - a*b + a*b*f0)a0 =a0a1 =a1a0 =a*b - a*b*f1a1 =b*f0 - a*b - b + a*b*f1a2 =a*b*f1 - a*b + 1a3 =b*f0 - b + 1a4 =b + 2*a*b - b*f0 - 2*a*b*f1 + 1A1 =[ b - b*f0, b - b*f1][ a*b - a - a*b*f0, a*b - a*b*f1]P22 =[ b, b*(b - b*f0) + a*b*(b - b*f1)] [ a*b, a*b*(a*b - a*b*f1) - b*(a - a*b + a*b*f0)] ans =2------闭环系统能控，能控性不变------A =0.5654 0.5654-1.4630 1.9034B =0.56541.9034C =1.0e+03 *0.27892.2767D =1ans =1.2344 + 0.6162i1.2344 - 0.6162ians =1.37961.3796------开环系统不稳定------ P2 =0.5654 1.39581.90342.7956 ans =2------开环系统能控------a0 =4285991457983477/2251799813685248 - (4285991457983477*f1)/2251799813685248a1 =(2827*f0)/5000 + (4285991457983477*f1)/2251799813685248 - 3474474420400697637/1407374883553280000ans =闭环系统稳定条件: 1-a0>0,1+a1+a0>0 ,1-a1+a0>0a2 =(4285991457983477*f1)/2251799813685248 - 2034191644298229/2251799813685248a3 =(2827*f0)/5000 + 2173/5000a4 =3780296982596825381/703687441776640000 - (4285991457983477*f1)/1125899906842624 - (2827*f0)/5000%打印稳定区域hold ongrid onlim=[-2,8,-1,3];fplot('[0.4746,1.4113-0.1485*x]',lim)plot([-0.7687 -0.7687],ylim,'m')shg运行结果：%验证闭环系统稳定性F=[1 0.5]L*eye(size(A))-(A-B*F)roots(poly(A-B*F))abs(roots(poly(A-B*F)))if abs(roots(poly(A-B*F)))<1disp('---------闭环系统稳定----------') else disp('----------闭环系统不稳定-----------') end运行结果F =1.0000 0.5000ans =[ lambda, -2827/10000] [ 2104/625, lambda - 4285991457983477/4503599627370496] ans =0.4758 + 0.8516i0.4758 - 0.8516ians =0.97550.9755---------闭环系统稳定----------结构图：六、设计使国民收入平稳增长的调控政策，并说明政策含义。