matlab求解动态规划
真的要开始动态规划了
基础部分
先看一下符号函数如何求值
clc;clear
syms f x1 x2
f=exp(-exp(-(x1+x2))) - x2*(1+x1^2);
symvar(f) %该函数返回的是符号函数中的自变量
g=matlabFunction(f);
g(1,1)
看下面这样的一个吃蛋糕问题
这里感觉原先教材上的写法很繁琐,不如用非线性规划来做;
有确切的状态转移方程,那非线性规划其实就是可以处理这种问题
这里用线性规划求解以下这一题,附上中间引用的函数
function f=cake(x)
beta=0.2;
f=-(log(x(1))+beta*log(x(2))+beta^2*log(x(3)));
end
% w=33%w0是蛋糕的分量
%
% %这里先考虑三天的情形,用非线性规划求解
% x0=[0;0;0]
% [x,fval]=fmincon(@cake,x0,[],[],[1 1 1],w,[0 0 0],[w w w]);
% x
% fval=-fval
satisfy=[];
for i=1:100x0=[0;0;0][x,fval]=fmincon(@cake,x0,[],[],[1 1 1],i,[0 0 0],[i i i]);xfval=-fvalsatisfy=[satisfy;fval];
end
plot(1:100,satisfy,'ob')%这里勾画出蛋糕的量与所获得最大满意度的关系图
得出的蛋糕的量和最优策略获得的满意度的变化值如上图所示
dynprog
接下来介绍一下dynprog函数的用法,这个函数是个人开发的一个函数,需要放进matlab toolbox 下面的datafun编辑器里面
具体的m文件可以在这个网址下载:http://www.verysource.com/code/30133134_1/DYNPROG.M.html
源码我顺带也直接附在下面了,需要的小伙伴可以自己存个m文件
function [p_opt,fval]=dynprog(x,DecisFun,ObjFun,TransFun)
% [p_opt,fval]=dynprog(x,DecisFun,ObjFun,TransFun)
% 自由始端和终端的动态规划,求指标函数最小值的逆序算法递归
% 计算程序。x是状态变量,一列代表一个阶段状态;M-函数
% DecisFun(k,x)由阶段k的状态变量x求出相应的允许决策变量;
% M-函数ObjFun(k,x,u)是阶段指标函数,M-函数TransFun(k,x,u)
% 是状态转移函数,其中x是阶段k的某状态变量,u是相应的决策变量;
% 输出p_opt由4列构成,p_opt=[序号组;最优策略组;最优轨线组;
% 指标函数值组];fval是一个列向量,各元素分别表示p_opt各
% 最优策略组对应始端状态x的最优函数值;
%
%例(参看胡良剑等编《数学实验--使用MATLAB》P180
%先写3个函数
% eg13f1_2.m
% function u=DecisF_1(k,x)
% 在阶段k由状态变量x的值求出其相应的决策变量所有的取值
% c=[70,72,80,76];q=10*[6,7,12,6];
% if q(k)-x<0,u=0:100; %决策变量不能取为负值
% else,u=q(k)-x:100;end; %产量满足需求且不超过100
% u=u(:);
% eg13f2_2.m
% function v=ObjF_1(k,x,u)
% 阶段k的指标函数
% c=[70,72,80,76];v=c(k)*u+2*x;
% eg13f3_2.m
% function y=TransF_1(k,x,u)
% 状态转移方程
% q=10*[6,7,12,6];y=x+u-q(k);
%调用DynProg.m计算如下:
% clear;x=nan*ones(14,4);% x是10的倍数,最大范围0≤x≤130,
% %因此x=0,1,...13,所以x初始化取14行,nan表示无意义元素
% x(1:7,1)=10*(0:6)'; % 按月定义x的可能取值
% x(1:11,2)=10*(0:10)';x(1:12,3)=10*(2:13)';
% x(1:7,4)=10*(0:6)';
% [p,f]=dynprog(x,'eg13f1_2','eg13f2_2','eg13f3_2')
% By X.D. Ding June 2000
k=length(x(1,:));f_opt=nan*ones(size(x));d_opt=f_opt;
t_vubm=inf*ones(size(x));x_isnan=~isnan(x);t_vub=inf;
% 计算终端相关值
tmp1=find(x_isnan(:,k));tmp2=length(tmp1);
for i=1:tmp2u=feval(DecisFun,k,x(i,k));tmp3=length(u);for j=1:tmp3tmp=feval(ObjFun,k,x(tmp1(i),k),u(j));if tmp<=t_vub, f_opt(i,k)=tmp;d_opt(i,k)=u(j);t_vub=tmp;
end;end;end
% 逆推计算各阶段的递归调用程序
for ii=k-1:-1:1tmp10=find(x_isnan(:,ii));tmp20=length(tmp10);for i=1:tmp20u=feval(DecisFun,ii,x(i,ii));tmp30=length(u);for j=1:tmp30tmp00=feval(ObjFun,ii,x(tmp10(i),ii),u(j));tmp40=feval(TransFun,ii,x(tmp10(i),ii),u(j));tmp50=x(:,ii+1)-tmp40;tmp60=find(tmp50==0);if ~isempty(tmp60),tmp00=tmp00+f_opt(tmp60(1),ii+1); if tmp00<=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(tmp1,1);
% 记录最优决策、最优轨线和相应指标函数值
p_opt=[];tmpx=[];tmpd=[];tmpf=[];
tmp0=find(x_isnan(:,1));tmp01=length(tmp0);
for i=1:tmp01,tmpd(i)=d_opt(tmp0(i),1); tmpx(i)=x(tmp0(i),1);tmpf(i)=feval(ObjFun,1,tmpx(i),tmpd(i));p_opt(k*(i-1)+1,[1,2,3,4])=[1,tmpx(i),...
tmpd(i),tmpf(i)];for ii=2:ktmpx(i)=feval(TransFun,ii-1,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(ObjFun,ii,tmpx(i),tmpd(i));p_opt(k*(i-1)+ii,[1,2,3,4])=[ii,tmpx(i),...
tmpd(i),tmpf(i)];
end;end;
下面用这道例题使用dynprog函数进行求解
首先要先准备三个函数
第一个decision函数,这里是对于每一次决策量的约束
第二个准备指标函数obj
第三个准备状态转移方程
接下来,再准备一下x(每一个阶段的初始状态)的矩阵,列对应这第n期,行对应着不同的可能取值
clc;clear
x=nan*ones(14,4);
x(1:7,1)=10*(0:6)';
x(1:11,2)=10*(0:10)';
x(1:12,3)=10*(2:13)';
x(1:7,4)=10*(0:6)';
看一下这里的起始矩阵长什么亚子
下面调用dynprog函数求解这个动态规划过程
[p,f]=dynprog(x,'decision','obj','trans')
p所对应的四列分别是 期数 存储量 生产策略 当季花销总额
p.s. 不同的期数 1 2 3 4 对应着不同的起始存储量
matlab求解动态规划相关推荐
- MATLAB数学建模-规划模型总结| MATLAB求解
目录 1 线性规划问题(LP) 风格1 风格2 2 非线性规划 3 动态规划 A星算法 基于dijkstra的概率路线图 4 多目标规划 帕累托最优 支配(Dominace) 不可支配解集 帕累托最优 ...
- MATLAB实现动态规划算法,基于Matlab的动态规划算法的实现及应用
陈甜甜 [摘要]介绍了动态规划的基本理论,包括动态规划的基本概念和基本原理,并针对生产与存储问题进行了分析,然后结合Matlab做了编程处理,使复杂问题简单化,从而使问题能更方便地得到解决. [关键词 ...
- matlab解符号方程组,matlab 求解符号方程组
1特殊符号可爱组成的小狗图案 求解符号方程组: 特殊符号可爱组成的小狗图案,缺失:matlab求解符号方程组4057/9 ▄██████▄ █████████▄ ███ ▄████▄▄▄▄███ ██ ...
- matlab微分方程组边值,matlab求解常微分方程边值问题的方法
matlab求解常微分方程边值问题的方法 Matlab 求解常微分方程边值问题的方法:bvp4c 函数 常微分方程的边值问题,即 boundary value problems ,简称 BVP 问题, ...
- 【数字信号处理】线性常系数差分方程 ( 使用 matlab 求解 “ 线性常系数差分方程 “ 示例 | A 向量分析 | B 向量分析 | 输入序列分析 | matlab 代码 )
文章目录 一.使用 matlab 求解 " 线性常系数差分方程 " 示例 1.B 向量元素 : x(n) 参数 2.A 向量元素 : y(n) 参数 3.输入序列 4.matlab ...
- 【数字信号处理】线性常系数差分方程 ( 卷积 与 “ 线性常系数差分方程 “ | 使用 matlab 求解 “ 线性常系数差分方程 “ )
文章目录 一.卷积 与 " 线性常系数差分方程 " 二.使用 matlab 求解 " 线性常系数差分方程 " 一.卷积 与 " 线性常系数差分方程 & ...
- 2021-01-13 Matlab求解微分代数方程 (DAE)
Matlab求解微分代数方程 (DAE) 什么是微分代数方程? 微分代数方程是一类微分方程,其中一个或多个因变量导数未出现在方程中.方程中出现的未包含其导数的变量称为代数变量,代数变量的存在意味着不能 ...
- MATLAB求解常微分方程
MATLAB求解微分方程_Falcon的博客-CSDN博客_matlab微分方程求解 matlab求解常微分方程(组)---dsolve.ode系列函数详解(含例程)_假电工的真的博客-CSDN博客_ ...
- 常微分方程matlab求解
常微分方程matlab求解 一般格式 matlab求解常微分方程的调用格式为: 例如,现在需要求解常微分方程 则有 y=dsolve('Dy=-2*y+2*x^2+2*x','x') 这个常微分方程的 ...
最新文章
- DllMain中不当操作导致死锁问题的分析--加载卸载DLL与DllMain死锁的关系
- Spring.NET性能
- 【Android 高性能音频】OboeTest 音频性能测试应用 ( 应用简介 | 测试内容 | 输出测试 | Oboe 缓冲区 与 工作负载修改 | 测试案例 )
- Asp.Net站点整合Discuz论坛实现同步注册和单点登录
- 【每日SQL打卡】​​​​​​​​​​​​​​​DAY 14丨重新格式化部门表【难度中等】
- angular 居中_Angular 的模块间通信
- 人工智能的未来-揭示人类思维的奥秘How to create a mind - Ray Kurzweil
- 复杂网络分析软件NetworkX和UCINET数据关联的方法
- 摄氏度符号英文计算机语言,温度表示-摄氏度怎样用英文表示温度?给几个例子,好吗? 爱问知识人...
- 计算机平面设计专业可以考什么证,平面设计师资格证怎么考_计算机平面设计职称...
- Codeforce Gym 100015I Identity Checker 暴力
- java bigdecimal 取整数_java-检查BigDecimal是否为整数值
- 用手机访问计算机共享资源,怎么进入共享文件夹?手机访问电脑局域网共享文件夹的方法...
- 2021年数学建模国赛湖北赛区推荐国奖名单
- java报表是什么_什么是报表工具
- Must call super constructor in derived class before accessing or returning from derived const
- java购物车设计_Java面向对象课程设计——购物车
- 巴菲特50年投资之道
- Java核心技术36讲(个人整理)
- Android 自定义View实现文本水平方向的跑马灯效果