如何做下列关于最优化的题目？

请查找数值分析

最小二乘法

在我们研究两个变量(x, y)之间的相互关系时，通常可以得到一系列成对的数据(x1, y1、x2, y2... xm , ym)；将这些数据描绘在x -y直角坐标系中(如图1), 若发现这些点在一条直线附近，可以令这条直线方程如(式1-1)。

Y计= a0 + a1 X (式1-1)

其中：a0、a1 是任意实数

为建立这直线方程就要确定a0和a1，应用《最小二乘法原理》，将实测值Yi与利用(式1-1)计算值(Y计=a0+a1X)的离差(Yi-Y计)的平方和〔∑(Yi - Y计)2〕最小为“优化判据”。

令: φ = ∑(Yi - Y计)2 (式1-2)

把(式1-1)代入(式1-2)中得:

φ = ∑(Yi - a0 - a1 Xi)2 (式1-3)

当∑(Yi-Y计)平方最小时，可用函数 φ 对a0、a1求偏导数，令这两个偏导数等于零。

(式1-4)

(式1-5)

亦即：

m a0 + (∑Xi ) a1 = ∑Yi (式1-6)

(∑Xi ) a0 + (∑Xi2 ) a1 = ∑(Xi, Yi) (式1-7)

得到的两个关于a0、 a1为未知数的两个方程组，解这两个方程组得出：

a0 = (∑Yi) / m - a1(∑Xi) / m (式1-8)

a1 = [∑Xi Yi - (∑Xi ∑Yi)/ m] / [∑Xi2 - (∑Xi)2 / m)] (式1-9)

这时把a0、a1代入(式1-1)中, 此时的(式1-1)就是我们回归的元线性方程即：数学模型。

在回归过程中，回归的关联式是不可能全部通过每个回归数据点(x1, y1、 x2, y2...xm,ym),为了判断关联式的好坏,可借助相关系数“R”，统计量“F”，剩余标准偏差“S”进行判断；“R”越趋近于 1 越好；“F”的绝对值越大越好；“S”越趋近于 0 越好。

R = [∑XiYi - m (∑Xi / m)(∑Yi / m)]/ SQR{[∑Xi2 - m (∑Xi / m)2][∑Yi2 - m (∑Yi / m)2]} (式1-10) ＊

在(式1-1)中，m为样本容量，即实验次数；Xi、Yi分别任意一组实验X、Y的数值。微积分应用课题一 最小二乘法

从前面的学习中, 我们知道最小二乘法可以用来处理一组数据, 可以从一组测定的数据中寻求变量之间的依赖关系, 这种函数关系称为经验公式. 本课题将介绍最小二乘法的精确定义及如何寻求 与 之间近似成线性关系时的经验公式. 假定实验测得变量之间的 个数据 , , …, , 则在 平面上, 可以得到 个点 , 这种图形称为“散点图”, 从图中可以粗略看出这些点大致散落在某直线近旁, 我们认为 与 之间近似为一线性函数, 下面介绍求解步骤.

考虑函数 , 其中 和 是待定常数. 如果 在一直线上, 可以认为变量之间的关系为 . 但一般说来, 这些点不可能在同一直线上. 记 , 它反映了用直线 来描述 , 时, 计算值 与实际值 产生的偏差. 当然要求偏差越小越好, 但由于 可正可负, 因此不能认为总偏差 时, 函数 就很好地反映了变量之间的关系, 因为此时每个偏差的绝对值可能很大. 为了改进这一缺陷, 就考虑用 来代替 . 但是由于绝对值不易作解析运算, 因此, 进一步用 来度量总偏差. 因偏差的平方和最小可以保证每个偏差都不会很大. 于是问题归结为确定 中的常数 和 , 使 为最小. 用这种方法确定系数 , 的方法称为最小二乘法.

由极值原理得 , 即

解此联立方程得

(*)

问题 I 为研究某一化学反应过程中, 温度 ℃)对产品得率 (％)的影响, 测得数据如下:

温度 ℃)

100 110 120 130 140 150 160 170 180 190

得率 (％)

45 51 54 61 66 70 74 78 85 89

(1) 利用“ListPlot”函数, 绘出数据 的散点图(采用格式: ListPlot[{ , , …, }, Prolog->AbsolutePointSize[3]] )

(2) 利用“Line”函数, 将散点连接起来, 注意观察有何特征? (采用格式: Show[Graphics[Line[{ , , …, }]] , Axes->True ])

(3) 根据公式(*), 利用“Apply”函数及集合的有关运算编写一个小的程序, 求经验公式

(程序编写思路为: 任意给定两个集合A (此处表示温度)、B(此处表示得率), 由公式(*)可定义两个二元函数(集合A和B为其变量)分别表示 和 . 集合A元素求和: Apply[Plus,A] 表示将加法施加到集合A上, 即各元素相加, 例如Apply[Plus,{1,2,3}]=6Length[A]表示集合A 元素的个数, 即为nA.B表示两集合元素相乘相加A*B表示集合A与B元素对应相乘得到的新的集合.)

(4) 在同一张图中显示直线 及散点图

(5) 估计温度为200时产品得率.

然而, 不少实际问题的观测数据 , , …, 的散点图明显地不能用线性关系来描叙, 但确实散落在某一曲线近旁, 这时可以根据散点图的轮廓和实际经验, 选一条曲线来近似表达 与 的相互关系.

问题 II 下表是美国旧轿车价格的调查资料, 今以 表示轿车的使用年数, (美元)表示相应的平均价格, 求 与 之间的关系.

使用年数

1 2 3 4 5 6 7 8 9 10

平均价格

2651 1943 1494 1087 765 538 484 290 226 204

(1) 利用“ListPlot”函数绘出数据 的散点图, 注意观察有何特征?

(2) 令 , 绘出数据 的散点图, 注意观察有何特征?

(3) 利用“Line”函数, 将散点 连接起来, 说明有何特征?

(4) 利用最小二乘法, 求 与 之间的关系

(5) 求 与 之间的关系

(6) 在同一张图中显示散点图 及 关于 的图形.

思考与练习

1. 假设一组数据 : , , …, 变量之间近似成线性关系, 试利用集合的有关运算, 编写一简单程序: 对于任意给定的数据集合 , 通过求解极值原理所包含的方程组, 不需要给出 、 计算的表达式, 立即得到 、 的值, 并就本课题 I /(3)进行实验.

注: 利用Transpose函数可以得到数据A的第一个分量的集合, 命令格式为:

先求A的转置, 然后取第一行元素, 即为数据A的第一个分量集合, 例如

(A即为矩阵 )

= (数据A的第一个分量集合)

= (数据A的第二个分量集合)

B-C表示集合B与C对应元素相减所得的集合, 如 = .

2. 最小二乘法在数学上称为曲线拟合, 请使用拟合函数“Fit”重新计算 与 的值, 并与先前的结果作一比较.

注: Fit函数使用格式:

设变量为x, 对数据A进行线性拟合, 如对题1中的A拟合函数为:

例1 求 f = 2 在0<x<8中的最小值与最大值

主程序为wliti1.m:

f='2*exp(-x).*sin(x)'

fplot(f,[0,8])%作图语句

[xmin,ymin]=fminbnd (f, 0,8)

f1='-2*exp(-x).*sin(x)'

[xmax,ymax]=fminbnd (f1, 0,8)

运行结果：

xmin = 3.9270ymin = -0.0279

xmax = 0.7854 ymax = 0.6448

★（借助课件说明过程、作函数的图形）

例2 对边长为3米的正方形铁板，在四个角剪去相等的正方形以制成方形无盖水槽，问如何剪法使水槽的容积最大？

设剪去的正方形的边长为x，则水槽的容积为： ，建立无约束优化模型为：min y=- ， 0<x<1.5

先编写M文件fun0.m如下:

function f=fun0(x)

f=-(3-2*x).^2*x

主程序为wliti2.m:

[x,fval]=fminbnd('fun0',0,1.5)

xmax=x

fmax=-fval

运算结果为: xmax = 0.5000,fmax =2.0000.即剪掉的正方形的边长为0.5米时水槽的容积最大,最大容积为2立方米.

★（借助课件说明过程、作函数的图形、并编制计算程序）

例3

1、编写M-文件 fun1.m:

function f = fun1 (x)

f = exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1)

2、输入M文件wliti3.m如下:

x0 = [-1, 1]

x=fminunc(‘fun1’,x0)

y=fun1(x)

3、运行结果:

x= 0.5000 -1.0000

y = 1.3029e-10

★（借助课件说明过程、作函数的图形并编制计算程序）

例4 Rosenbrock 函数 f（x1，x2）=100（x2-x12）2+(1-x1)2 的最优解（极小）为x*=（1，1），极小值为f*=0.试用不同算法（搜索方向和步长搜索）求数值最优解.初值选为x0=（-1.2 , 2）.

为获得直观认识，先画出Rosenbrock 函数的三维图形, 输入以下命令：

[x,y]=meshgrid(-2:0.1:2,-1:0.1:3)

z=100*(y-x.^2).^2+(1-x).^2

mesh(x,y,z)

画出Rosenbrock 函数的等高线图,输入命令：

contour(x,y,z,20)

hold on

plot(-1.2,2,' o ')

text(-1.2,2,'start point')

plot(1,1,'o')

text(1,1,'solution')

f='100*(x(2)-x(1)^2)^2+(1-x(1))^2'

[x,fval,exitflag,output]=fminsearch(f, [-1.2 2])

运行结果:

x =1.00001.0000

fval =1.9151e-010

exitflag = 1

output =

iterations: 108

funcCount: 202

algorithm: 'Nelder-Mead simplex direct search'

★（借助课件说明过程、作函数的图形并编制计算程序）

（五）、 作业

陈酒出售的最佳时机问题

某酒厂有批新酿的好酒,如果现在就出售,可得总收入R0=50万元(人民币),如果窖藏起来待来日(第n年)按陈酒价格出售,第n年末可得总收入 (万元),而银行利率为r=0.05,试分析这批好酒窖藏多少年后出售可使总收入的现值最大. (假设现有资金X万元,将其存入银行,到第n年时增值为R(n)万元,则称X为R(n)的现值.)并填下表：

最数下降法解无约束优化的程序

先建立一维搜索的m文件：minWP.m如下

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x,minf]=minWP(f,XMAX,c1,c2,alpha,tol)

%一维搜索的Wolfe-Powell法

%作者:龚纯 王正林<<精通 Matlab 最优化计算>>

%f:目标函数

%XMAX:搜索最大值

%c1:可接受系数1

%c2:可接受系数2

%alpha:增大步长倍数

%tol:精度

%x:极小值点

%minf:极小值点处的函数值

format long

if nargin==5

tol=1.0e-6

end

if ~(c1>0)||~(c1<c2)||~(c2<1)||~(XMAX>0)||~(alpha>1)

error('参数不对')

end

var=findsym(f)

df=diff(f)

f0=subs(f,var,0)

df0=subs(df,var,0)

a=0

b=XMAX

if b<inf

t=(a+b)/2

else

t=10

end

while 1

ft=subs(f,var,t)

fu=f0+c1*t*df0

if ft<=fu

dft=subs(df,var,t)

dfl=c2*t*df0

if dft>dfl

x=t

break

else

a=t

if b==inf

t=t*alpha

else

t=(a+b)/2

end

continue

end

else

b=t

t=(a+b)/2

continue

end

end

minf=subs(f,var,t)

format short

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

然后建立最速下降法的m文件如下：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%minFD.m

function [x,minf]=minFD(f,x0,var,tol)

%最速下降法

%作者：龚纯 王正林《精通Matlab最优化计算》

%目标函数：f

%初始点：x0

%自变量向量：var

%精度：tol

%所求的驻点：x

%驻点处的函数值

format long

if nargin==3

tol=1.0e-6

end

gradf=jacobian(f,var) %f的梯度

wucha=1

syms lamda

while wucha>tol

d=-subs(gradf,var,x0)

wucha=norm(d)

y=x0+lamda*d

yf=subs(f,var,y)

%l=minHJ(yf,0,100)

l=minWP(yf,10,0.2,0.6,2)

x1=x0+l*d

x0=x1

end

x=x1

minf=subs(f,var,x)

format short

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

保存后就可调用了：

>>clear

>>syms x1 x2

>>f=1.5*x1^2+0.5*x2^2-x1*x2-2*x1

f =

3/2*x1^2+1/2*x2^2-x1*x2-2*x1

>>[x,mf]=minFD(f,[-2 4],[x1,x2])

x =

1.00001.0000

mf =

-1.0000

---