数学建模与数学实验

数学建模与数学实验数学建模与数学实验 MATLAB 入门入门 MATLAB 作为线性系统的一种分析和仿真工具 是理工科大学生应该掌握的技术 工具 它作为一种编程语言和可视化工具 可解决工程 科学计算和数学学科中许多 问题 MATLAB 建立在向量 数组和矩阵的基础上 使用方便 人机界面直观 输出结 果可视化 矩阵是 MATLAB 的核心 MATLAB 的进入与运行方式 两种 一 变一 变 量量 与与 函函 数数 1 变量 变量 MATLAB 中变量的命名规则是 1 变量名必须是不含空格的单个词 2 变量名区分大小写 3 变量名最多不超过 19 个字符 4 变量名必须以字母打头 之后可以是任意字母 数字或下划线 变量名中不 允许使用标点符号 特殊变量表特殊变量表 特殊 变量 取 值 ans用于结果的缺省变量名 pi圆周率 eps计算机的最小数 当和 1 相加就产生一个比 1 大的数 flops浮点运算数 inf无穷大 如 1 0 NaN不定量 如 0 0 i ji j 1 nargin所用函数的输入变量数目 nargou t 所用函数的输出变量数目 realmi n 最小可用正实数 realma x 最大可用正实数 2 数学运算符号及标点符号 数学运算符号及标点符号 加法运算 适用于两个数或两个同阶矩阵相 加 减法运算 乘法运算 点乘运算 除法运算 点除运算 乘幂运算 点乘幂运算 反斜杠表示左除 1 MATLAB 的每条命令后 若为逗号或无标点符号 则显示命令的结果 若命令 后为分号 则禁止显示结果 2 后面所有文字为注释 3 表示续行 3 数学函数 数学函数 函 数名 称函 数名 称 sin x 正弦函数asin x 反正弦函数 cos x 余弦函数acos x 反余弦函数 tan x 正切函数atan x 反正切函数 abs x 绝对值max x 最大值 min x 最小值sum x 元素的总和 sqrt x 开平方exp x 以 e 为底的指 数 log x 自然对数 log10 x 以 10 为底的对 数 sign x 符号函数fix x 取整 4 M 文件文件 MATLAB 的内部函数是有限的 有时为了研究某一个函数的各种性态 需要为 MATLAB 定义新函数 为此必须编写函数文件 函数文件是文件名后缀为 M 的文件 这类文件的第一行必须是一特殊字符 function 开始 格式为 function 因变量名 函数名 自变量名 函数值的获得必须通过具体的运算实现 并赋给因变量 M 文件建立方法 1 在 Matlab 中 点 File New M file 2 在编辑窗口中输入程序内容 3 点 File Save 存盘 M 文件名必须与函数名一致 Matlab 的应用程序也以 M 文件保存 例 例 定义函数 f x1 x2 100 x2 x12 2 1 x1 2 1 建立 M 文件 fun m function f fun x f 100 x 2 x 1 2 2 1 x 1 2 2 可以直接使用函数 fun m 例如 计算 f 1 2 只需在 Matlab 命令窗口键入命令 x 1 2 fun x 二 数二 数 组组 1 创建简单的数组 创建简单的数组 x a b c d e f 创建包含指定元素的行向量 x first last 创建从 first 开始 加 1 计数 到 last 结束的行向量 x first increment last 创建从 first 开始 加 increment 计数 last 结束的行向 量 x linspace first last n 创建从 first 开始 到 last 结束 有 n 个元素的行向量 x logspace first last n 创建从开始 到结束 有 n 个元素的对数分隔行向量 2 数组元素的访问数组元素的访问 1 访问一个元素访问一个元素 x i 表示访问数组 x 的第 i 个元素 2 访问一块元素访问一块元素 x a b c 表示访问数组 x 的从第 a 个元素开始 以步长为 b 到第 c 个元素 但不超过 c b 可以为负数 b 缺损时为 1 3 直接使用元素编址序号直接使用元素编址序号 x a b c d 表示提取数组 x 的第 a b c d 个元 素构成一个新的数组 x a x b x c x d 3 数组的方向 数组的方向 前面例子中的数组都是一行数列 是行方向分布的 称之为行向量 数组也可以 是列向量 它的数组操作和运算与行向量是一样的 唯一的区别是结果以列形式显示 产生列向量有两种方法 直接产生 例 c 1 2 3 4 转置产生 例 b 1 2 3 4 c b 说明 以空格或逗号分隔的元素指定的是不同列的元素 而以分号分隔的元素指定了 不同行的元素 4 数组的运算 数组的运算 1 标量 标量 数组运算数组运算 数组对标量的加 减 乘 除 乘方是数组的每个元素对该标量施加相应的加 减 乘 除 乘方运算 设 设 a a1 a2 an c 标量标量 则 则 a c a1 c a2 c an c a c a1 c a2 c an c a c a1 c a2 c an c 右除 右除 a c c a1 c a2 c an 左除 左除 a c a1 c a2 c an c c a c a1 c a2 c an 2 数组 数组 数组运算数组运算 当两个数组有相同维数时 加 减 乘 除 幂运算可按元素对元素方式进行的 不同大小或维数的数组是不能进行运算的 设 设 a a1 a2 an b b1 b2 bn 则 则 a b a1 b1 a2 b2 an bn a b a1 b1 a2 b2 an bn a b a1 b1 a2 b2 an bn a b b1 a1 b2 a2 bn an a b a1 b1 a2 b2 an bn 三 三 矩矩 阵阵 1 矩阵的建立 矩阵的建立 逗号或空格用于分隔某一行的元素 分号用于区分不同的行 除了分号 在输入矩 阵时 按 Enter 键也表示开始一新行 输入矩阵时 严格要求所有行有相同的列 例 m 1 2 3 4 5 6 7 8 9 10 11 12 p 1 1 1 1 2 2 2 2 3 3 3 3 特殊矩阵的建立特殊矩阵的建立 a 产生一个空矩阵 当对一项操作无结果时 返回空矩阵 空矩阵的大小为零 b zeros m n 产生一个 m 行 n 列的零矩阵 c ones m n 产生一个 m 行 n 列的元素全为 1 的矩阵 d eye m n 产生一个 m 行 n 列的单位矩阵 2 矩阵中元素的操作 矩阵中元素的操作 1 矩阵 A 的第 r 行 A r 2 矩阵 A 的第 r 列 A r 3 依次提取矩阵 A 的每一列 将 A 拉伸为一个列向量 A 4 取矩阵 A 的第 i1 i2 行 第 j1 j2 列构成新矩阵 A i1 i2 j1 j2 5 以逆序提取矩阵 A 的第 i1 i2 行 构成新矩阵 A i2 1 i1 6 以逆序提取矩阵 A 的第 j1 j2 列 构成新矩阵 A j2 1 j1 7 删除 A 的第 i1 i2 行 构成新矩阵 A i1 i2 8 删除 A 的第 j1 j2 列 构成新矩阵 A j1 j2 9 将矩阵 A 和 B 拼接成新矩阵 A B A B 3 矩阵的运算 矩阵的运算 1 标量 矩阵运算 同标量 数组运算 2 矩阵 矩阵运算 1 元素对元素的运算 同数组 数组运算 2 矩阵运算 矩阵加法 A B 矩阵乘法 A B 方阵的行列式 det A 方阵的逆 inv A 方阵的特征值与特征向量 V D eig A 关系与逻辑运算关系与逻辑运算 1 关系操作符 关系操作符 关系操关系操 作符作符 说明说明 小于 大于 大于或等于 等于 不等于 2 逻辑运算符逻辑运算符 逻辑操作符逻辑操作符说明说明 与 或 非 控制流控制流 MATLAB 提供三种决策或控制流结构 for 循环 while 循环 if else end 结构 这些结构经常包含大量的 MATLAB 命令 故经常出现在 MATLAB 程序中 而不是直接加 在 MATLAB 提示符下 1 for 循环 允许一组命令以固定的和预定的次数重复循环 允许一组命令以固定的和预定的次数重复 for x array commands end 在 for 和 end 语句之间的命令串 commands 按数组 array 中的每一列执行一次 在 每一次迭代中 x 被指定为数组的下一列 即在第 n 次循环中 x array n 例例 对 n 1 2 10 求 xn sin npi 10 的值 2 While 循环循环 与与 for 循环以固定次数求一组命令相反 循环以固定次数求一组命令相反 while 循环以不定的次数求一组语句的值循环以不定的次数求一组语句的值 while expression commands end 只要在表达式 expression 里的所有元素为真 就执行 while 和 end 语句之间的命令串 commands 例例 设银行年利率为 11 25 将 10000 元钱存入银行 问多长时间会连本带利翻一番 3 If Else End 结构结构 1 有一个选择的一般形式是 有一个选择的一般形式是 if expression commands end 如果在表达式 expression 里的所有元素为真 就执行 if 和 end 语句之间的命令串 commands 1 2 12 11 2 ff xx xx xf求设例 先建立 M 文件 fun1 m 定义函数 f x 再在 Matlab 命令窗口输入 fun1 2 fun1 1 即可 2 有三个或更多的选择的一般形式是 有三个或更多的选择的一般形式是 if expression1 commands1 else if expression2 commands2 else if expression3 commands3 else if else commands end end end end 1 5 0 2 0 x 1x02 1x1x f x 3 2 fff x x求设例 先建立 M 文件 fun2 m 定义函数 f x 再在 Matlab 命令窗口输入 fun2 2 fun2 0 5 fun2 1 即 可 实验作业实验作业 对以下问题 编写 M 文件 1 用起泡法对 10 个数由小到大排序 即将相邻两个数比较 将小的调到前头 2 有一个矩阵 编程求出其最大值及其所处的位置 3 编程求一球从 100 米高度自由落下 每次落地后反跳回原高度的一半 再落下 求它在第 10 次落地时 共经过多少米 第 10 次反弹有多高 4 有一函数 写一程序 输入自变量的值 输出函数值 2sin 2 xyxyxf MATLAB 作图作图 二二 维维 图图 形形 1 曲线图曲线图 Matlab 作图是通过描点 连线来实现的 故在画一个曲线图形之前 必须先取得该图形上 的一系列的点的坐标 即横坐标和纵坐标 然后将该点集的坐标传给 Matlab 函数画图 命令为 PLOT X Y S X Y 是向量是向量 分别表示点集的横坐标和纵坐标分别表示点集的横坐标和纵坐标 S 是线型是线型 PLOT X Y 画实线画实线 PLOT X Y1 S1 X Y2 S2 X Yn Sn 将多条线画在将多条线画在 一起一起 例例 在 0 2 pi 用红线画 sin x 用绿圈画 cos x 解解 x linspace 0 2 pi 30 y sin x z cos x plot x y r x z g0 2 符号函数符号函数 显函数 隐函数和参数方程显函数 隐函数和参数方程 画图画图 1 ezplot ezplot f x a b 表示在 a x b 绘制显函数 f f x 的函数图 ezplot f x y xmin xmax ymin ymax 表示在区间 xmin x xmax 和 ymin y ymax 绘 制 隐函数 f x y 0 的函数图 ezplot x t y t tmin tmax 表示在区间 tmin t tmax 绘制参数方程 x x t y y t 的函数图 例例 在 0 pi 上画 y cos x 的图形 解解 输入命令 ezplot sin x 0 pi 例例 在 0 2 pi 上画 星形图tx 3 cos ty 3 sin 解解 输入命令 ezplot cos t 3 sin t 3 0 2 pi 例例 在 2 0 5 0 2 上画隐函数的图0 sin xyex 解解 输入命令 ezplot exp x sin x y 2 0 5 0 2 2 fplot fplot fun lims 表示绘制字符串 fun 指定的函数在 lims xmin xmax 的图形 注意 1 fun 必须是 M 文件的函数名或是独立变量为 x 的字符串 2 fplot 函数不能画参数方程和隐函数图形 但在一个图上可以画多个图形 例例 在 1 2 上画的 图形 3sin 22 xey x 解解 先建 M 文件 myfun1 m function Y myfun1 x Y exp 2 x sin 3 x 2 再输入命令 fplot myfun1 1 2 例例 在 2 2 范围内绘制函数 tanh 的图形 解解 fplot tanh 2 2 例例 x y 的取值范围都在 2 2 画函数 tanh x sin x cos x 的图形 解解 输入命令 fplot tanh x sin x cos x 2 pi 1 1 1 1 3 对数坐标图对数坐标图 在很多工程问题中 通过对数据进行对数转换可以更清晰地看出数据的某些特征 在对数 坐标系中描绘数据点的曲线 可以直接地表现对数转换 对数转换有双对数坐标转换和单轴对 数坐标转换两种 用 loglog 函数可以实现双对数坐标转换 用 semilogx 和 semilogy 函数可以 实现单轴对数坐标转换 loglog Y 表示 x y 坐标都是对数坐标系 semilogx Y 表示 x 坐标轴是对数坐标系 semilogy 表示 y 坐标轴是对数坐标系 plotyy 有两个 y 坐标轴 一个在左边 一个在右边 例例 用方形标记创建一个简单的 loglog 解解 输入命令 x logspace 1 2 loglog x exp x s grid on 标注格栅 例例 创建一个简单的半对数坐标图 解解 输入命令 x 0 1 10 semilogy x 10 x 例例 绘制 y x3 的函数图 对数坐标图 半对数坐标图 三三 维维 图图 形形 空空 间间 曲曲 线线 1 一条曲线一条曲线 PLOT3 x y z s x y z n 维向量 分别表示曲线上点集的横坐标 纵坐标 函数值 s 指定颜色 线形等 例例 在区间 0 10 pi 画出参数曲线 x sin t y cos t z t 解解 t 0 pi 50 10 pi plot3 sin t cos t t rotate3d 旋转 2 多条曲线 多条曲线 PLOT3 x y z 其中 x y z 是都是 m n 矩阵 其对应的每一列表示一条曲线 例 画多条曲线观察函数 Z X Y 2 解 x 3 0 1 3 y 1 0 1 5 X Y meshgrid x y Z X Y 2 plot3 X Y Z 这里 meshgrid x y 的作用是产生一个以向量 x 为行 向量 y 为列的矩阵 空空 间间 曲曲 面面 1 surf x y z 画出数据点 x y z 表示的曲面 数据矩阵 x y z 分别表示数据点的横坐标 纵坐标 函数值 例 画函数 Z X Y 2 的图形 解 x 3 0 1 3 y 1 0 1 5 X Y meshgrid x y Z X Y 2 surf X Y Z shading flat 将当前图形变得平滑 2 Mesh x y z 画网格曲面 数据矩阵 分别表示数据点的横坐标 纵坐标 函数值 例 画出曲面 Z X Y 2 在不同视角的网格图 解 x 3 0 1 3 y 1 0 1 5 X Y meshgrid x y Z X Y 2 mesh X Y Z 3 meshz X Y Z 在网格周围画一个 curtain 图 如 参考平面 例 绘 peaks 的网格图 解 输入命令 X Y meshgrid 3 125 3 Z praks X Y Meshz X Y Z 处处 理理 图图 形形 1 在图形上加格栅 图例和标注 在图形上加格栅 图例和标注 1 GRID ON 加格栅在当前图上 GRID OFF 删除格栅 2 hh xlabel string 在当前图形的 x 轴上加图例 string hh ylabel string 在当前图形的 y 轴上加图例 string hh zlabel string 在当前图形的 z 轴上加图例 string hh title string 在当前图形的顶端上加图例 string 例例 在区间 0 2 pi 画 sin x 的图形 并加注图例 自变量 X 函数 Y 示意图 并加格 栅 解解 x linspace 0 2 pi 30 y sin x plot x y xlabel 自变量 X ylabel 函数 Y title 示意图 grid on 3 hh gtext string 命令 gtext string 用鼠标放置标注在现有的图上 运行命令 gtext string 时 屏幕上 出现当前图形 在图形上出现一个交叉的十字 该十字随鼠标的移动移动 当按下鼠标左 键时 该标注 string 放在当前十交叉的位置 例例 在区间 0 2 pi 画 sin x 并分别标注 sin x cos x 解解 x linspace 0 2 pi 30 y sin x z cos x plot x y x z gtext sin x gtext cos x 2 定制坐标 定制坐标 Axis xmin xmax ymin ymax zmin zmax 定制图形坐标 x y z 的最大 最小值 Axis auto 将坐标轴返回到自动缺省值 例例 在区间 0 005 0 01 显示 sin 1 x 的图形 解解 x linspace 0 0001 0 01 1000 y sin 1 x plot x y axis 0 005 0 01 1 1 3 图形保持 图形保持 1 hold on 保持当前图形 以便继续画图到当前图上 hold of 释放当前图形窗口 例例 将 y sin x y cos x 分别用点和线画出在同一屏幕上 解解 x linspace 0 2 pi 30 y sin x z cos x plot x z hold on Plot x y 2 figure h 新建 h 窗口 激活图形使其可见 并把它置于其它图形之上 例例 区间 0 2 pi 新建两个窗口分别画出 y sin x z cos x 解解 x linspace 0 2 pi 100 y sin x z cos x plot x y title sin x pause figure 2 plot x z title cos x 4 分割窗口 分割窗口 h subplot mrows ncols thisplot 划分整个作图区域为 mrows ncols 块 逐行对块访问 并激活第 thisplot 块 其后的作 图语句将图形画在该块上 subplot mrows ncols thisplot 激活已划分为 mrows ncols 块的屏幕中的第 thisplot 块 其后的作图语句将图形画在块 上 subplot 1 1 1 命令 Subplot 1 1 1 返回非分割状态 例例 将屏幕分割为四块 并分别画出 y sin x z cos x a sin x cos x b sin x cos x 解解 x linspace 0 2 pi 100 y sin x z cos x a sin x cos x b sin x cos x eps subplot 2 2 1 plot x y title sin x subplot 2 2 2 plot x z title cos x subplot 2 2 3 plot x a title sin x cos x subplot 2 2 4 plot x b title sin x cos x 5 缩放图形 缩放图形 zoom on 为当前图形打开缩放模式 单击鼠标左键 则在当前图形窗口中 以鼠标点中的点为中心的图形放大 2 倍 单击鼠 标右键 则缩小 2 倍 zoom off 关闭缩放模式 例例 缩放 y sin x 的图形 解解 x linspace 0 2 pi 30 y sin x Plot x y zoom on 6 改变视角改变视角 view 1 view a b 命令 view a b 改变视角到 a b a 是方位角 b 为仰角 缺省视角为 37 5 30 2 view x y z view 用空间矢量表示的 三个量只关心它们的比例 与数值的大小无关 x 轴 view 1 0 0 y 轴 view 0 1 0 z 轴 view 0 0 1 例 画出曲面 Z X Y 2 在不同视角的网格图 解 x 3 0 1 3 y 1 0 1 5 X Y meshgrid x y Z X Y 2 subplot 2 2 1 mesh X Y Z subplot 2 2 2 mesh X Y Z view 50 34 subplot 2 2 3 mesh X Y Z view 60 70 subplot 2 2 4 mesh X Y Z view 0 1 1 7 动画动画 Moviein getframe movie 函数 Moviein 产生一个帧矩阵来存放动画中的帧 函数 getframe 对当前的图象进行快 照 函数 movie 按顺序回放各帧 例 将曲面 peaks 做成动画 解 x y z peaks 30 surf x y z axis 3 3 3 3 10 10 m moviein 15 for i 1 15 view 37 5 24 i 1 30 m i getframe end movie m 特特 殊殊 的的 二二 维维 图图 形形 函函 数数 1 极坐标图 极坐标图 polar theta rho s 用角度 theta 弧度表示 和极半径 rho 作极坐标图 用 s 指定线型 例 的极坐标图形 2cos2sin r 解 theta linspace 0 2 pi rho sin 2 theta cos 2 theta polar theta rho g title Polar plot of sin 2 theta cos 2 theta 2 散点图散点图 scatter X Y S C 在向量 X 和 Y 的指定位置显示彩色圈 X 和 Y 必须大小相同 例例 绘制 seamount 散点图 解解 输入命令 load seamount scatter x y 5 z 3 平面等值线图 平面等值线图 contour x y z n 绘制 n 个等值线的二维等值线图 例例 在范围 2 x 2 2 y 3 内绘的等值线图 22 yx xez 解解 输入命令 X Y meshgeid 2 2 2 2 2 3 Z X exp X 2 Y 2 C h contour X Y Z clabel C h colormap cool 特特 殊殊 的的 三三 维维 图图 形形 函函 数数 1 空间等值线图 空间等值线图 contour 3 x y z n 其中 n 表示等值线数 例例 山峰的三维和二维等值线图 解 x y z peaks subplot 1 2 1 contour3 x y z 16 s grid xlabel x axis ylabel y axis zlabel z axis title contour3 of peaks subplot 1 2 2 contour x y z 16 s grid xlabel x axis ylabel y axis title contour of peaks 绘制山区地貌图绘制山区地貌图 要在某山区方圆大约 27 平方公里范围内修建一条公路 从山脚出发经过一个居民区 再到达一个矿区 横向纵向分别每隔 400 米测量一次 得到一些地点的高程 平面区域 0 x 5600 0 y0 使B km 1 中至少有一个大于零 km 1 km 2 kmm 2 所有的 0 1 2 m i i 则一定存在另一个可行基 它对应的基可行解使目标函数值更小 改进方法改进方法 令 则把从原有的基中取出来 把加进后得到的 0 kmi i kmi 0 min kml l l x km x 仍是基 即是所要找的新基 1 x 2 x l x km x 1 l x m x 用用 MATLAB 优化工具箱解线性规划优化工具箱解线性规划 1 模型模型 min z cX bAXts 命令命令 x linprog c A b 2 模型 模型 min z cX bAXts beqXAeq 命令命令 x linprog c A b Aeq beq 注意 若没有不等式 存在 则令 A b 3 模型模型 min z cX VLB X VUBbAXts beqXAeq 命令命令 1 x linprog c A b Aeq beq VLB VUB 2 x linprog c A b Aeq beq VLB VUB X0 注意 1 若没有等式约束 则令 Aeq beq beqXAeq 2 其中 X0 表示初始点 4 命令 x fval linprog 返回最优解 及 处的目标函数值 fval 例例 1 max 654321 6 064 0 72 0 32 0 28 0 4 0 xxxxxxz 85003 0 03 0 03 0 01 0 01 0 01 0 654321 xxxxxxts 70005 0 02 0 41 xx 10005 0 02 0 52 xx 90008 0 03 0 63 xx6 2 10 jxj 解解 编写 M 文件 xxgh1 m 如下 c 0 4 0 28 0 32 0 72 0 64 0 6 A 0 01 0 01 0 01 0 03 0 03 0 03 0 02 0 0 0 05 0 0 0 0 02 0 0 0 05 0 0 0 0 03 0 0 0 08 b 850 700 100 900 Aeq beq vlb 0 0 0 0 0 0 vub x fval linprog c A b Aeq beq vlb vub 例例 2 321 436minxxxz 120 321 xxxts 30 1 x 500 2 x 20 3 x 3 2 1 436 min x x x z 50 120 010 111 3 2 1 x x x ts 3 2 1 20 0 30 x x x 解解 编写 M 文件 xxgh2 m 如下 c 6 3 4 A 0 1 0 b 50 Aeq 1 1 1 beq 120 vlb 30 0 20 vub x fval linprog c A b Aeq beq vlb vub 例例 3 问题一的解答 改写为 Xz8121110913min S t 900 800 3 12 15 0000 00011 14 0 X 500 600 400 100100 010010 001001 X 0 6 5 4 3 2 1 x x x x x x X 编写 M 文件 xxgh3 m 如下 f 13 9 10 11 12 8 A 0 4 1 1 1 0 0 0 0 0 0 0 5 1 2 1 3 b 800 900 Aeq 1 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 1 beq 400 600 500 vlb zeros 6 1 vub x fval linprog f A b Aeq beq vlb vub 结果 x 0 0000 600 0000 0 0000 400 0000 0 0000 500 0000 fval 1 3800e 004 即在甲机床上加工 600 个工件 2 在乙机床上加工 400 个工件 1 500 个工件 3 可在满足条 件的情况下使总加工费最小为 13800 例例 2 问题二的解答 改写为 s t 2 1 3640min x x z 45 35 2 1 x x 编写 M 文件 xxgh4 m 如下 c 40 36 A 5 3 b 45 Aeq beq vlb zeros 2 1 vub 9 15 调用 linprog 函数 x fval linprog c A b Aeq beq vlb vub 结果为 x 9 0000 0 0000 fval 360 即只需聘用 9 个一级检验员 注 注 本问题应还有一个约束条件 x1 x2 取整数 故它是一个整数线性规划问题 这里把 它当成一个线性规划来解 求得其最优解刚好是整数 x1 9 x2 0 故它就是该整数规划 的最优解 若用线性规划解法求得的最优解不是整数 将其取整后不一定是相应整数规划 的最优解 这样的整数规划应用专门的方法求解 投资的收益和风险投资的收益和风险 一 问题提出一 问题提出 市场上有 n 种资产 i 1 2 n 可以选择 现用数额为 M 的相当大的资金作一个 i s 时期的投资 这 n 种资产在这一时期内购买的平均收益率为 风险损失率为 投 i s i r i q 资越分散 总的风险越小 总体风险可用投资的中最大的一个风险来度量 i s 购买时要付交易费 费率 当购买额不超过给定值时 交易费按购买计算 i s i p i u i u 另外 假定同期银行存款利率是 既无交易费又无风险 5 0 r 0 r 已知 n 4 时相关数据如下 i s i r i q i p 元 i u S1282 51103 S2211 52198 S3235 54 552 S4252 66 540 试给该公司设计一种投资组合方案 即用给定达到资金 M 有选择地购买若干种资产或存 银行生息 使净收益尽可能大 使总体风险尽可能小 二 基本假设和符号规定二 基本假设和符号规定 基本假设 基本假设 1 投资数额 M 相当大 为了便于计算 假设 M 1 2 投资越分散 总的风险越小 3 总体风险用投资项目中最大的一个风险来度量 i s 4 n 种资产 S 之间是相互独立的 i 5 在投资的这一时期内 ri pi qi r0为定值 不受意外因素影响 6 净收益和总体风险只受 ri pi qi影响 不受其他因素干扰 符号规定 符号规定 Si 第 i 种投资项目 如股票 债券 ri pi qi 分别为 Si的平均收益率 风险损失率 交易费率 ui Si的交易定额 同期银行利率 0 r xi 投资项目 Si的资金 a 投资风险度 Q 总体收益 Q 总体收益的增量 三 模型的建立与分析三 模型的建立与分析 1 总体风险用所投资的 Si 中最大的一个风险来衡量 即 max qixi i 1 2 n 2 购买 Si所付交易费是一个分段函数 即 pixi xi ui 交易费 piui xi ui 而题目所给定的定值 ui 单位 元 相对总投资 M 很小 piui更小 可以忽略不计 这样购买 Si的净收益为 ri pi xi 3 要使净收益尽可能大 总体风险尽可能小 这是一个多目标规划模型 目标函数 MAX n i iii xpr 0 MINmax qixi 约束条件 M n i ii xp 0 1 xi 0 i 0 1 n 4 模型简化 a 在实际投资中 投资者承受风险的程度不一样 若给定风险一个界限 a 使最大的一个 风险 qixi M a 可找到相应的投资方案 这样把多目标规划变成一个目标的线性规划 模型模型 1 1 固定风险水平 优化收益 目标函数 Q MAX 1 1 n i iii xpr 约束条件 a M xq ii xi 0 i 0 1 n Mxp ii 1 b 若投资者希望总盈利至少达到水平 k 以上 在风险最小的情况下寻找相应的投资组合 模型模型 2 固定盈利水平 极小化风险 目标函数 R min max qixi 约束条件 k n i iii xpr 0 xi 0 i 0 1 n Mxp ii 1 投资者在权衡资产风险和预期收益两方面时 希望选择一个令自己满意的投资组合 因此对风险 收益赋予权重 s 0 s 1 s 称为投资偏好系数 模型模型 3 目标函数 min s max qixi 1 s n i iii xpr 0 约束条件 M xi 0 i 0 1 2 n n i ii xp 0 1 四 模型四 模型 1 的求解的求解 模型 1 为 minf 0 05 0 27 0 19 0 185 0 185 x0 x1 x2 x3 x 4 T x0 1 01x1 1 02x2 1 045x3 1 065x4 1 s t 0 025x1 a 0 015x2 a 0 055x3 a 0 026x4 a xi 0 i 0 1 4 由于 a 是任意给定的风险度 到底怎样给定没有一个准则 不同的投资者有不同的风险度 我们从 a 0 开始 以步长 a 0 001 进行循环搜索 编制程序如下 a 0 while 1 1 a 1 c 0 05 0 27 0 19 0 185 0 185 Aeq 1 1 01 1 02 1 045 1 065 beq 1 A 0 0 025 0 0 0 0 0 0 015 0 0 0 0 0 0 055 0 0 0 0 0 0 026 b a a a a vlb 0 0 0 0 0 vub x val linprog c A b Aeq beq vlb vub a x x Q val plot a Q axis 0 0 1 0 0 5 hold on a a 0 001 end xlabel a ylabel Q 计算结果 a 0 0030 x 0 4949 0 1200 0 2000 0 0545 0 1154 Q 0 1266 a 0 0060 x 0 0 2400 0 4000 0 1091 0 2212 Q 0 2019 a 0 0080 x 0 0000 0 3200 0 5333 0 1271 0 0000 Q 0 2112 a 0 0100 x 0 0 4000 0 5843 0 0 Q 0 2190 a 0 0200 x 0 0 8000 0 1882 0 0 Q 0 2518 a 0 0400 x 0 0000 0 9901 0 0000 0 0 Q 0 2673 五 五 结果分析结果分析 1 风险大 收益也大 2 当投资越分散时 投资者承担的风险越小 这与题意一致 即 冒险的投资者会出现集中 投资的情况 保守的投资者则尽量分散投资 3 曲线上的任一点都表示该风险水平的最大可能收益和该收益要求的最小风险 对于不同 风险的承受能力 选择该风险水平下的最优投资组合 4 在 a 0 006 附近有一个转折点 在这一点左边 风险增加很少时 利润增长 很快 在这一点右边 风险增加很大时 利润增长很缓慢 所以对于风险和 收益没有特殊偏好的投资者来说 应该选择曲线的拐点作为最优投资组合 大约是 a 0 6 Q 20 所对应投资方案为 风险度 收益 x0 x1 x2 x3 x4 0 0060 0 2019 0 0 2400 0 4000 0 1091 0 2212 实验作业实验作业 某厂生产甲乙两种口味的饮料 每百箱甲饮料需用原料 6 千克 工人 10 名 可获利 10 万元 每 百箱乙饮料需用原料 5 千克 工人 20 名 可获利 9 万元 今工厂共有原料 60 千克 工人 150 名 又由于其他条件所限甲饮料产量不超过 8 百箱 问如何安排生产计划 即两种饮料各生产多少 使获利最大 进一步讨论 1 若投资 0 8 万元可增加原料 1 千克 问应否作这项投资 2 若每百箱甲饮料获利可增加 1 万元 问应否改变生产计划 无约束最优化无约束最优化 求解无约束最优化问题的基本思想求解无约束最优化问题的基本思想 标准形式标准形式 其中 Xf n EX min 1 EEf n Xf n EX max minXf n EX 求解的基本思想求解的基本思想 以二元函数为例 连续可微 2 x 21x xf 0 1 x 2 x 0 5 3 1 0 X 1 X 2 X 0 Xf 1 Xf 2 Xf 搜索过程搜索过程 最优点 1 1 初始点 1 1 2 1 22 1221 1 100 minxxxxxf 1 x 2 x f 114 00 0 790 583 39 0 530 232 60 0 180 001 50 0 09 0 030 98 0 370 110 47 0 590 330 20 0 800 630 05 0 950 900 003 0 990 991E 4 0 9990 9981E 5 0 99970 99981E 8 无约束优化问题的基本算法无约束优化问题的基本算法 1 最速下降法 共轭梯度法 算法步骤 最速下降法 共轭梯度法 算法步骤 给定初始点 允许误差 令 k 0 n EX 0 0 计算 k Xf 检验是否满足收敛性的判别准则 k Xf 若满足 则停止迭代 得点 否则进行 k XX 令 从出发 沿进行一维搜索 kk XfS k X k S 即求使得 k k k kkk SXfSXf 0 min 令 k k 1 返回 k k kk SXX 1 最速下降法是一种最基本的算法 它在最优化方法中占有重要地位 最速下降法的优点是工 作量小 存储变量较少 初始点要求不高 缺点是收敛慢 最速下降法适用于寻优过程的前 期迭代或作为间插步骤 当接近极值点时 宜选用别种收敛快的算法 2 牛顿算法步骤牛顿算法步骤 1 选定初始点 给定允许误差 令 k 0 n EX 0 0 2 求 检验 若 则停止迭代 否则转向 3 k Xf 1 2 k Xf k Xf k XX 3 令 牛顿方向 kkk XfXfS 12 4 转回 2 kkk SXX 1 1 kk 如果 f 是对称正定矩阵 A 的二次函数 则用牛顿法经过一次迭代 就可达到最优点 如不是 二次函数 则牛顿法不能一步达到极值点 但由于这种函数在极值点附近和二次函数很近 似 因此牛顿法的收敛速度还是很快的 牛顿法的收敛速度虽然较快 但要求 Hessian 矩阵要可逆 要计算二阶导数和逆矩阵 就 加大了计算机计算量和存储量 3 拟牛顿法 拟牛顿法 为克服牛顿法的缺点 同时保持较快收敛速度的优点 利用第 k 步和第 k 1 步得到的 k X 构造一个正定矩阵近似代替 或用近 1 k X k Xf 1 k Xf 1 k G 2k Xf 1 k H 似代替 将牛顿方向改为 12 k Xf 1 k G 1 k S 1 k Xf 1 k S 1 k H 1 k Xf 从而得到下降方向 通常采用迭代法计算 迭代公式迭代公式为 1 k G 1 k H BFGSBFGS Boryden Fletcher Goldfarb Shanno 公式 kkTk kTkkk kTk Tkk kk xGx GxxG xf ff GG 1 kTk Tkk kTk kkTk kk xf xx xf fHf HH 1 1 kTk TkkkkTkk xf xfHHfx DFPDFP Davidon Fletcher Powell 公式 kTk Tkk kTk kkTk kk Xf ff fX XGX GG 1 1 kTk TkkkkTkk fX fXGGXf 计算时可置 单位阵 对于给出 kkTk kTkkk kTk Tkk kk fHf HffH Xf XX HH 1 IH 1 的利用上面的公式进行递推 这种方法称为拟牛顿法拟牛顿法 1 X Matlab 优化工具箱简介优化工具箱简介 1 MATLAB 求解优化问题的主要函数求解优化问题的主要函数 类 型模 型基本函数名 一元函数极小Min F x s t x1 x x2 x fminbnd F x1 x2 无约束极小 Min F X X fminunc F X0 X fminsearch F X0 线性规划 Min XcT s t AX b X linprog c A b 二次规划 Min xTHx cTx 2 1 s t Ax b X quadprog H c A b 约束极小 非线性规划 Min F X s t G X 0 X fmincon FG X0 达到目标问题 Min r s t F x wr goal X fgoalattain F x goal w 极小极大问题 Min max Fi x X Fi x s t G x 0 则x 为解 否则 x不是最终解 它只是迭代制止 时优化过程的值 所有优化函数 fval 解x处的目标函数值 linprog quadprog fgoalattain fmincon fminimax lsqcurvefit lsqnonlin fminbnd exitflag 描述退出条件 exitflag 0 表目标函数收敛于解x处 exitflag 0 表已达到函数评价或迭代 的最大次数 exitflag 0 表目标函数不收敛 output 包含优化结果信息的输出结构 所有优化函数 Iterations 迭代次数 Algorithm 所采用的算法 FuncCount 函数评价次数 控制参数控制参数 options 可以通过函数可以通过函数 optimset 创建或修改 命令的格式如下 创建或修改 命令的格式如下 1 options optimset optimfun 创建一个含有所有参数名 并与优化函数 optimfun 相关的默认值的选项结构 options 2 options optimset param1 value1 param2 value2 创建一个名称为 options 的优化选项参数 其中指定的参数具有指定值 所有未指定的参数 取默认值 3 options optimset oldops param1 value1 param2 value2 创建名称为 oldops 的参数的拷贝 用指定的参数值修改 oldops 中相应的参数 例 opts optimset Display iter TolFun 1e 8 该语句创建一个称为 opts 的优化选项结构 其中显示参数设为 iter TolFun 参数设为 1e 8 用用 Matlab 解无约束优化问题解无约束优化问题 1 一元函数无约束优化问题一元函数无约束优化问题 min f x 21 xxx 常用格式如下 常用格式如下 1 x fminbnd fun x1 x2 2 x fminbnd fun x1 x2 options 3 x fval fminbnd 4 x fval exitflag fminbnd 5 x fval exitflag output fminbn