基于Matlab的动态规划程序实现

最优控制实验2大连海事大学自动化专业动态规划方法的Matlab实现与应用动态规划(DynamicProgramming)是求解决策过程最优化的有效数学方法，它是根据“最优决策的任何截断仍是最优的”这最优性原理，通过将多阶段决策过程转化为一系列单段决策问题，然后从最后一段状态开始逆向递推到初始状态为止的一套最优化求解方法。1．动态规划基本组成(1)阶段整个问题的解决可分为若干个阶段依次进行，描述阶段的变量称为阶段变量，记为k(2)状态状态表示每个阶段开始所处的自然状况或客观条件，它描述了研究问题过程的状况。各阶段状态通常用状态变量描述，用kx表示第k阶段状态变量，n个阶段决策过程有n+1个状态。(3)决策从一确定的状态作出各种选择从而演变到下一阶段某一状态，这种选择手段称为决策。描述决策的变量称为决策变量，决策变量限制的取值范围称为允许决策集合。用()kkux表示第k阶段处于状态kx时的决策变量，它是kx的函数。用()kkDxDk(xk)表示kx的允许决策的集合。(4)策略每个阶段的决策按顺序组成的集合称为策略。由第k阶段的状态kx开始到终止状态的后部子过程的策略记为{}11(),(),,()kkkknnuxuxux++。可供选择的策略的范围称为允许策略集合，允许策略集合中达到最优效果的策略称为最优策略。从初始状态*11()xx=出发，过程按照最优策略和状态转移方程演变所经历的状态序列{}****121,,,,nnxxxx+称为最优轨线。(5)状态转移方程如果第k个阶段状态变量为kx，作出的决策为ku，那么第k+1阶段的状态变量1kx+也被完全确定。用状态转移方程表示这种演变规律，记为1(,)kkkxTxu+=。(6)指标函数指标函数是系统执行某一策略所产生结果的数量表示，是衡量策略优劣的数量指标，它定义在全过程和所有后部子过程上，用()kkfx表示。过程在某阶段j的阶段指标函数是衡量该阶段决策优劣数量指标，取决于状态jx和决策ju，用(,)jjjvxu表示。2．动态规划基本方程(){}11()min,,(),()kkkkkkkkkkfxgvxufxuDx++=∈⎡⎤⎣⎦Matlab实现(dynprog.m文件)function[p_opt,fval]=dynprog(x,DecisFun,SubObjFun,TransFun,ObjFun)%x是状态变量，一列代表一个阶段的所有状态；%M-函数DecisFun(k,x)由阶段k的状态变量x求出相应的允许决策变量;%M-函数SubObjFun(k,x,u)是阶段指标函数,%M-函数ObjFun(v,f)是第k阶段至最后阶段的总指标函数%M-函数TransFun(k,x,u)是状态转移函数,其中x是阶段k的某状态变量,u是相应的决策变量;%输出p_opt由4列构成，p_opt=[序号组;最优策略组;最优轨线组;指标函数值组]；%输出fval是一个列向量，各元素分别表示p_opt各最优策略组对应始端状态x的最优函数值。最优控制实验2大连海事大学自动化专业k=length(x(1,:));%判断决策级数x_isnan=~isnan(x);%非空状态矩阵t_vubm=inf*ones(size(x));%性能指标中间矩阵f_opt=nan*ones(size(x));%总性能指标矩阵d_opt=f_opt;%每步决策矩阵tmp1=find(x_isnan(:,k));%最后一步状态向量tmp2=length(tmp1);%最后一步状态个数fori=1:tmp2u=feval(DecisFun,k,x(tmp1(i),k));tmp3=length(u);%决策变量forj=1:tmp3%求出当前状态下所有决策的最小性能指标tmp=feval(SubObjFun,k,x(tmp1(i),k),u(j));iftmp=t_vubm(i,k)%t_vubf_opt(i,k)=tmp;d_opt(i,k)=u(j);t_vubm(i,k)=tmp;end;end;endforii=k-1:-1:1tmp10=find(x_isnan(:,ii));tmp20=length(tmp10);fori=1:tmp20%求出当前状态下所有可能的决策u=feval(DecisFun,ii,x(tmp10(i),ii));tmp30=length(u);forj=1:tmp30%求出当前状态下所有决策的最小性能指标tmp00=feval(SubObjFun,ii,x(tmp10(i),ii),u(j));%单步性能指标tmp40=feval(TransFun,ii,x(tmp10(i),ii),u(j));%下一状态tmp50=x(:,ii+1)-tmp40;%找出下一状态在x矩阵的位置tmp60=find(tmp50==0);if~isempty(tmp60),ifnargin5,tmp00=tmp00+f_opt(tmp60(1),ii+1);%setthedefaultobjectvalueelse,tmp00=feval(ObjFun,tmp00,f_opt(tmp60(1),ii+1));end%当前状态的性能指标iftmp00=t_vubm(i,ii)f_opt(i,ii)=tmp00;d_opt(i,ii)=u(j);t_vubm(i,ii)=tmp00;end;end;end;end;end;fval=f_opt(:,1);tmp0=find(~isnan(fval));fval=fval(tmp0,1);p_opt=[];tmpx=[];tmpd=[];tmpf=[];tmp01=length(tmp0);fori=1:tmp01tmpd(i)=d_opt(tmp0(i),1);tmpx(i)=x(tmp0(i),1);tmpf(i)=feval(SubObjFun,1,tmpx(i),tmpd(i));p_opt(k*(i-1)+1,[1,2,3,4])=[1,tmpx(i),tmpd(i),tmpf(i)];forii=2:ktmpx(i)=feval(TransFun,ii,tmpx(i),tmpd(i));tmp1=x(:,ii)-tmpx(i);tmp2=find(tmp1==0);if~isempty(tmp2),tmpd(i)=d_opt(tmp2(1),ii);end;tmpf(i)=feval(SubObjFun,ii,tmpx(i),tmpd(i));p_opt(k*(i-1)+ii,[1,2,3,4])=[ii,tmpx(i),tmpd(i),tmpf(i)];end;end;最优控制实验2大连海事大学自动化专业3．动态规划方法应用例某电子设备由5种元件1,2,3,4,5组成，其可靠性分别为0.9,0.8,0.5,0.7,0.6为保证电子设备系统的可靠性，同种元件可并联多个。现允许设备使用元件的总数为15个，问如何设计使设备可靠性最大。将该问题看成一个5阶段动态规划问题，每个元件的配置看成一个阶段。记kx为配置第k个元件时可用元件的总数(状态变量)；ku为第k个元件并联的数目(决策变量);kc为第k个元件的可靠性，阶段指标函数为(,)1(1)kukkkkvxuc=−−;状态转移方程为1kkkxxu+=−基本方程为(){}{}111()min,,()minkkkkkkkkkfxgvxufxvf+++==⎡⎤⎣⎦i根据上述的允许决策、阶段指标函数、状态转移方程和基本方程写出下面的4个M-函数%DecisF1.mfunctionu=DecisF1(k,x)%在阶段k由状态变量x的值求出其相应的决策变量所有的取值ifk==5,u=x;elseu=1:x-1;end%SubObjF1.mfunctionv=SubObjF1(k,x,u)%阶段k的指标函数c=[0.9,0.8,0.5,0.5,0.4];v=1-(1-c(k)).^u;v=-v;%将求max转换为求min%TransF1.mfunctiony=TransF1(k,x,u)%状态转移方程y=x-u;%ObjF1.mfunctiony=ObjF1(v,f)%基本方程中的函数gy=v*f;y=-y;%将求max转换为求min调用DynProg.m计算的主程序如下(example1.m):clear;n=15;%15个元件x1=[n;nan*ones(n-1,1)];x2=1:n;x2=x2’;x=[x1,x2,x2,x2,x2]；[p,f]=dynprog(x,’DecisF1’,’SubObjF1’,’TransF1’,’ObjF1’)运行结果：p=1.000015.00002.0000-0.99002.000013.00002.0000-0.96003.000011.00004.0000-0.93754.00007.00003.0000-0.97305.00004.00004.0000-0.9744f=-0.8447结果表明1,2,3,4和5号元件分别并联2,2,4,3,4个，系统总可靠性最大为0.8447最优控制实验2大连海事大学自动化专业4．实验要求(1)根据前面2，3部分的说明，编写出dynprog.m，DecisF1.m，SubObjF1.m，TransF1.m，ObjF1.m五个文件，并运行example1.m文件，查看结果是否正确。(2)现有4种不同的车床1,2,3,4,同时加工500件相同的零件。各车床加工一个零件的时间分别为0.5,0.1,0.2和0.05.问如何给4个车床分配加工零件数目,使完工时间最短?要求参考(1)编写出DecisF2.m，SubObjF2.m，TransF2.m，ObjF2.m四个文件，通过调用dynprog函数，求解本问题。[提示：可将问题按车床编号分为4个阶段。设状态变量kx表示分配给第k号车床至第4号车床的零件总数。决策变量ku表示分配给第k号车床的零件数,1,2,,kkux=且有44ux=.状态转移方程1kkkxxu+=−,阶段指标函数()kkvu表示ku个零件分配到第k号车床加工所需的时间()kkkkvuut=,kt是k号车床的加工时间.()kkfx表示kx个零件分配给第k至第4号车床加工所需的最短时间(用时最长的车床所需时间为总加工时间)](3)已知离散系统方程x(k+1)=x(k)+u(k),x(1)=1,}1,1{)(−∈ku,4321()()()()()kJxkukxkukxk=⎡⎤=++⎣⎦∑,求控制序列u*(k)和状态x*(k),使J=min.要求参考(1)编写出DecisF3.m，SubObjF3.m，TransF3.m，ObjF3.m四个文件，通过调用dynprog函数，求解本问题。