学的时候没太在意，不过在后面遇到了蛮多相关的知识点（机器学习、起重机等都有它的身影）。

概述

优化的概念

可以理解为：求极小值
优化问题的解决离不开数学手段
优化技术要解决的问题
- 数学建模
- 优化问题的求解
- 优化结果的评价

最优化问题的数学模型

例 : 某工厂生产A、B产品如下表，合理分配使利润最大

产品	材料	工时	电力	利润
A	9	3	4	60
B	4	10	5	120
供应量	360	300	200	-

优化问题的数学模型的概念

设计变量：设计方案中待求的可变的量（是独立变量)，用 $X=(x_1,x_2…x_n)^T$ 表示
预定参数：设计中的已知项
设计空间：以设计分量 $x_i$ 为坐标轴所构成的空间
目标函数：将优化的目标表达成设计变量的函数，记为 $F(X)$，用于评价方案的好坏
约束函数限制条件：用 $g_u(X)\leq0,h_v(X)=0(u,v=1,2…)$
可行域：同时满足 $g_u(X),h_v(X)$ 的点的集合
等值线：等值线上的所有 $X$ 使目标函数或约束函数值相等

优化问题的数学模型的标准模式

求设计变量 $X=(x_1,x_2…x_n)^T$ ，使得

$f(X)\rightarrow\min$

并且满足

$g_u(X)\leq0,u=1,2...\\ h_v(X)=0,v=1,2...$

优化问题的分类

有无约束
目标函数约束函数性质
- 线性规划：$F,g,h$ 中均为线性
- 非线性规划：$F,g,h$ 中 有一个不是 线性
目标函数个数
- 单目标优化：模型中只有一个目标函数
- 多目标优化：模型有多个目标函数
设计变量

优化问题的求解思路

解析法
数值计算方法

对于优化问题，依次从初始点出发，找到一串点 $X^0,X^1…X^k$ ，对应的目标函数值 $f(X^0)>f(X^1)>…$ 若存在最优解，$f(X)$ 最终收敛于 $f(X^)$ ，$X^$ 为最优解

根据 $X^k$ 和 $X^*$ 的接近程度，有 3 种收敛条件

$|X^{k+1}-X^k|<\varepsilon$
$|F(X^{k+1})-F(X^k)|<\varepsilon$
$|\nabla F(X^k)|<\varepsilon$

迭代一般格式

给定初始点 $X^{k}|_{k=0}$，收敛精度 $\varepsilon$
按某种规则构造搜索方向 $S^{k}$
从 $X^{k}$ 出发 , 沿 $S^k$ 方向，找到在此方向上的最优点 $X^{k+1}$

x[k]={0};
for(k=0;Z(X)<ε;k++){
    X[k+1]=X[k]+α[k]S[k]
}

其中Z(X) 为收敛条件，α[k]为步长因子

数学基础

极值理论

一元函数极值

必要条件：$f’(x_0)=0$ ，$x_0$ 为驻点
充分条件：$f’’(x_0)>0$ ， $x_0$ 为极小值；$f’’(x_0)<0$ ， $x_0$ 为极大值；

二元函数极值

必要条件：$f’_x(x_0,y_0)=0,f’_y(x_0,y_0)=0$
充分条件：若 $\partial^2f/\partial x^2=A$，$\partial^2f/\partial xy=B$，$\partial^2f/\partial x^2=C$，$AC-B^2>0$ 有极值；$AC-B^2<0$ 无极值

多元函数极值

方向导数与梯度

以二元函数为例，在点 $X^0$ 沿任意方向 $S_0$ 的变化率

$\frac{\partial f}{\partial S_0}\\ =\lim\limits_{ds\rightarrow0}\frac{f(X^1)-f(X^0)}{ds}\\ =\lim\limits_{dx,dy\rightarrow0}\left[\frac{f(x_0+dx,y_0)-f(x_0,y_0)}{dx}\frac{dx}{ds}+\frac{f(x_0+dx,y_0+dy)-f(x_0+dx,y_0)}{dy}\frac{dy}{ds}\right]\\ =\frac{\partial f}{\partial x}\cos\theta_1+\frac{\partial f}{\partial y}\cos\theta_2=\begin{bmatrix} \frac{\partial f}{\partial x}\\ \frac{\partial f}{\partial y} \end{bmatrix}\begin{bmatrix} \cos\theta_1\\ \cos\theta_2 \end{bmatrix}=\nabla f\begin{bmatrix} \cos\theta_1\\ \cos\theta_2 \end{bmatrix}$

这里写成 $\cos\theta_1,\cos\theta_2$ 是考虑到多元函数

其中：$\nabla f$ 为梯度，$\begin{bmatrix}
\cos\theta_1\\
\cos\theta_2
\end{bmatrix}$ 为单位向量，记为 $S$

$\frac{\partial f}{\partial S}=|\nabla f|\cos\theta,\theta=<\nabla f,S>$

梯度性质

$\nabla F$ 表示该函数上升最快的方向，也是等值线的法线方向
梯度是函数在一点邻域的局部描述，离开邻域不一定是下降最快的方向

多元函泰勒展开

作用：对复杂函数在一点邻域内简化

$f(X)=f(X_0)+\nabla f(X)(X-X_0)+\frac12(X-X_0)^T\nabla^2f(X)(X-X_0)$

其中 $\nabla=\begin{bmatrix}
\partial/\partial x \\
\partial /\partial y
\end{bmatrix},\nabla^2=\begin{bmatrix}
\partial^2/\partial x^2 &\partial^2/\partial x\partial y\\
\partial ^2/\partial y\partial x&\partial^2/\partial y^2
\end{bmatrix}$

$\nabla^2f$ 又称海塞矩阵

函数在 $X^*$ 处取极值条件

必要条件：$\nabla f(X^*)=0$

对任意 $S,\nabla f·S\geq0$ ，取 $S=S_0,S=-S_0$ 得到
充分条件：$\nabla^2f(X^*)$ 正定

其他点都要比极小点大

设 $M$ 是 $n$ 阶方阵，如果对任何非零向量 $z$ ，都有 $z^TMz\geq0$ ，其中 $z^T$ 表示 $z$ 的转置，就称 $M$ 为正定矩阵。

【补充：判断正定矩阵的方法】

计算对称矩阵A的各阶顺序主子式。若A的各阶顺序主子式均大于零，则A是正定的

函数凹凸性

凸集：任意连接点集中任意两点 $x_1,x_2$ 的线段全包含在集合内
凸函数：$f[\alpha x_1+(1-\alpha)x_2]\leq\alpha f(x_1)+(1-\alpha)f(x_2),0<\alpha<1$
若 $f(X)$ 的海塞矩阵 $\nabla f$ 处处正定，$f(X)$ 为严格凸函数
凹凸性应用
- 约束函数 $g(X)$ 是凸函数，$g(X)<0$ 围成区域为凸集
- 约束函数 $g(X)$ 是凹函数，$g(X)>0$ 围成区域为凸集
- 凹凸性与极值关系
  - 目标函数是凸函数，约束函数是凸函数，极值点在可行域内（a）
  - 目标函数是凸函数，约束函数是凸函数，极值点不在可行域内（b）
  - 目标函数是凸函数，约束函数不是凸函数（c）
  - 目标函数不是凸函数，约束函数是凸函数（d）

一维搜索

对于一个优化问题 $\min f(X)$ ，通过搜索方向 $S$ 的构造，从初始点 $X_0$ 出发，寻找一组是目标函数下降的点的序列。在迭代的每一步 $X^{k+1}=X^K+\alpha_KS^K$ ，带入到目标函数，得 $f(\alpha_K)$ ，每一步迭代都归结于求 $\alpha_k$ （$X,S$ 已知）

$f(X^K)\rightarrow f(X^{K+1})$ 转变为 $f(\alpha_K)\rightarrow f(X^{K=1})$

求解基本步骤

找到一个包含极小点的区间
在区间内找到 $\alpha^*$

进退法（搜索区间的确定）

从某给定初值 $\alpha_0$ 开始，以初始步长 $h$ 向前试探。如果函数值上升，则步长变号，即改变试探方向；如果函数值下降，则维持原来的试探方向。区间的始点、中间点依次沿试探方向移动一步。此过程一直进行到函数值再次上升时为止，即可找到搜索区间的终点。最后得到的三点即为搜索区间的始点、中间点和终点，形成函数值的 “高—低—高” 趋势。

步骤：给定初始点及步长 $\alpha_0,h$

若 $f(\alpha_k+h)>f(\alpha_k)$ ，最优步长位于 $\alpha_k+h$ 左侧，右端点可以确定。搜索方向不对，使 $h=-h$ ,重新搜索，直到出现 $f(a_K+h)>f(\alpha_k))$ 。这个 $\alpha_K+h$ 就是左端点
1
2
3
4
5
6
7
if (f(alpha_k + h) > f(alpha_k)){ b = alpha_k + h; while(f(alpha_k - h) > f(alpha_k)){ alpha_k -= h; } a = alpha_k - h; }

反过来也是类似的

if (f(alpha_k + h) < f(alpha_k)){
    a = alpha_k;
    while(f(alpha_k + h) > f(alpha_k)){
        alpha_k += h;
    }
    a = alpha_k + h;
}

进退法

区间消去法（不断减少区间最终获得极小点）

已知区间 $[a,b]$ 在其中任取 $a_1<b_2$

若 $f(a_1)<f(b_1)$ ，极值点在 $b_1$ 左侧，消去无用区间得到$[a,b_1]$
若 $f(a_1)=f(b_1)$ ，极值点在 $a_1,b_1$ 中间，消去无用区间得到$[a_1,b_1]$
若 $f(a_1)>f(b_1)$ ，极值点在 $a_1$ 右侧，消去无用区间得到$[a_1,b]$

黄金分割法

以搜索区间的0.618倍作为比较点 $a_1,b_1$ 的确定原则

$a_1=a+0.618(b-a)$
$b_1=b-0.618(b-a)$

黄金分割法

格点法

其算法为：给定区间 $[a,b]$，将其等分为 $n+1$ 个区间 $(n\geq3)$，计算 $f(a_1),f(a_2)…f(a_n)$，找出其中的最小点 $f(a_m)=\min\{f(a_1),f(a_2)…f(a_n)\}$，$a_m$左右两点就构成搜索区间 $[a_{m-1},a_{m+1}]$，直至 $|a_{m+1}-a_{m-1}|\leq\varepsilon$

格点法

平分法

取 $a,b$ 中点 $a_K=(a+b)/2$ 利用中点附近的导数的正负性，判断取舍区间

平分法

当一阶导数无法求时，可以使用 $f(\alpha_k+\varepsilon)-f(\alpha_k-\varepsilon)$ 代替

二次插值法

原理：

对函数$f(x)$取三点 ${x_1}<{x_2}<{x_3}$ 对应函数值 $f(x_1)>f(x_2)<f(x_3) $ ，以此构造插值多项式 $ P(x)=a_0+a_1x+a_2x^2 $ ，其极值点 $ x^*=-a_1/2a_2 $

二次插值法

无约束规划

上一章确定了步长，但是没有确定方向，这一章就是在无约束问题中确定方向的

梯度法（最速下降法）

$-\nabla f(X)$ 方向是函数值下降最快的方向，取方向 $S=-\nabla f(X)$

迭代格式：$X^{K+1}=X^K-\alpha_K\nabla f(X)$
收敛条件：$||\nabla f(X)||<\varepsilon$

特点

相邻两次搜索正交，搜索路径呈锯齿形，越靠近极点，越慢，搜索效率越低

牛顿法

由于高次函数的导数不好求，于是现将其用泰勒公式转化为二次 ，再求梯度

迭代格式 $X^{K+1}=X^K-[\nabla^2 f(X)]^{-1}\nabla f(X)$

流程图太简单了。。。。

特点

工作量大、对高次函数，泰勒展开只是近似解，且步长恒为1，其极值不准确

修正牛顿法

在牛顿法的基础上引入了最优步长，确保沿 $\nabla f(X)$ 方向下降最快

迭代格式 $X^{K+1}=X^K-\alpha_K[\nabla^2 f(X)]^{-1}\nabla f(X)$

流程图与梯度法类似，知识把迭代格式改了一下

共轭方向法

构造共轭方向，函数值下降的方向就是共轭方向，这样得出极小点

基本概念

设 $A$ 为正定矩阵，若一组向量 $s_1,s_2…s_k$ 满足 $s_i^TAs_j=0,i,j=1,2…k,i\neq j$ 则称这组向量关于 $A$ 共轭。当 $A=E$ 时，则这组向量为共轭向量。

设 $M$ 是 $n$ 阶方阵，如果对任何非零向量 $z$ ，都有 $z^TMz\geq0$ ，其中 $z^T$ 表示 $z$ 的转置，就称 $M$ 为正定矩阵。

共轭性质

要使 $X^2$ 为二元函数的极小点，只需要使两次一维搜索方向 $S^0,S^1$ 与 $\nabla^2 f(X)$ 共轭
对 n 元二次函数，若 $S^i,i=1,2…n$ 与 $\nabla^2f(X)$ 共轭，沿着这几个方向进行一维搜索，最多进行 n 次即可找到极小点

共轭方向的构造

平行搜索法

对函数 $f(X)=\frac12X^T\nabla^2fX+BX+C$ 任意两个点 $X_1^0,X_2^0$ ，分别沿着同一方向 $S^0$ 进行一维搜索，得到两个一维搜索点 $X_1^1,X_2^1$ ，连接两点的方向 $S^1=X_2^1-X_1^1$ 与 $S^0$ 关于 $\nabla^2f(X)$ 共轭

共轭梯度法

$S^0=\nabla f(X^0),S^{K+1}=-\nabla f(X)+\beta_KS^K,\beta_K=||\nabla f(X^{K+1})||^2/||\nabla f(X^{K})||^2$

实验报告

1. 进退法求解初始区间

实验日期：2018 年 11 月 8 日

1.1 实验目的

加深对进退法的基本理论和算法步骤的理解。
培养学生独立编制计算机程序的能力。
掌握常用进退法程序的使用方法。

1.2 上机内容

1.2.1 算法概述

进退法用于求解一元函数 $f(\alpha)$ 极小点 ${\alpha^❀}$ 所在区间 $[a,b]$ 。假设函数 $f(\alpha)$ 为单态函数，根据单态区间的性质，为了确定极小点 $\alpha^❀$ 所在区间 $[a,b]$ ，应使函数 $f(\alpha)$ 在 $[a,b]$ 区间形成 “高—低—高” 趋势。

其算法为：从某给定初值 $\alpha_0$ 开始，以初始步长 $h$ 向前试探。如果函数值上升，则步长变号，即改变试探方向；如果函数值下降，则维持原来的试探方向。区间的始点、中间点依次沿试探方向移动一步。此过程一直进行到函数值再次上升时为止，即可找到搜索区间的终点。最后得到的三点即为搜索区间的始点、中间点和终点，形成函数值的 “高—低—高” 趋势。

1.2.2 程序框图

1.2.3 源程序

/*************************************************
* @Author: Fu_Qingchen
* @Date: 2018-11-7 
* @Description: 进退法求解区间
*************************************************/
#include<stdio.h>
#include<stdlib.h>
#define A0 1234
#define H 1069

/*************************************************
* @Description: 要找区间的函数
*************************************************/
double function(double x) {
	return (x * x + 15 * x + 1);
}

/*************************************************
* @Description: 进退法主函数
* @Input: a0:初始值,h:步长,*range:区间左右端点
* @Return: 区间左右端点
*************************************************/
void jintui(double a0, double h, double *range) {
	if (function(a0)==function(a0+h))
	{
		*(range) = a0;
		*(range + 1) = a0 + h;
	}
	else if (function(a0) > function(a0 + h))
	{
		*(range) = a0;
		do
		{
			a0 += h;
		} while (function(a0)>function(a0 + h));
		*(range + 1) = a0 + h;
	}
	else
	{
		*(range + 1) = a0 + h;
		do
		{
			a0 -= h;
		} while (function(a0)<function(a0 + h));
		*(range) = a0 - h;
	}
}

void main() {
	double a0 = A0, h = H;
	double range[2] = {0,0};
	jintui(a0, h, range);
	printf("Result:%f,%f\n", range[0], range[1]);
	system("pause");
}

1.3 过程分析

遇到的问题：

1.3.1 数组在函数中的传递

我发现C语言函数一般只能返回一个值，传递数组就比较麻烦

解决方法

使用指针在函数之间传递数组

例如：主函数中，调用的 jintui 函数通过 range 传递区间左右两个端点值

1 2	`double range[2] = {0,0}; jintui(a0, h, range);`

1.4 结果分析

给定函数：

$f(\alpha)=\alpha^2+15\alpha+1$

其理论极小点为 $\alpha^*=-7.5$

给定初值 $\alpha_0=1234$ ，给定步长 $h=1069$

结果：

所求区间 $[-1973.000000,2303.000000]$ 包含极小点 $\alpha=-7.5$ ，通过验证

2. 一维搜索法求解最优步长

实验日期：2018 年 11 月 15 日

2.1 实验目的

加深对一维搜索方法的基本理论和算法步骤的理解。
培养学生独立编制计算机程序的能力。
掌握常用一维搜索方法程序的使用方法。

2.2 上机内容

2.2.1 算法概述

此处使用格点法进行求解。格点法用于通过给定区间 $[a,b]$ 求解函数极小点 $\alpha^*$ 。

2.2.2 程序框图

2.2.3 源程序

/*************************************************
* @Author: Fu_Qingchen
* @Date: 2018-10-30
* @Description: 格点法求解迭代步长
*************************************************/
#include<stdio.h>
#include<stdlib.h>
#define A0 999
#define H 10	//格点法步长
#define N 10 //格点法等分数
#define V 0.00001	//精度

/*************************************************
* @Description: 要找极小点的函数
*************************************************/
double function(double x) {
	return (x * x + 2 * x + 1);
}


/*************************************************
* @Description: 进退法主函数
* @Input: a0:初始值,h:步长,*range:区间左右端点
* @Return: 区间左右端点
*************************************************/
void jintui(double a0, double h, double *range) {
	if (function(a0) == function(a0 + h))
	{
		*(range) = a0;
		*(range + 1) = a0 + h;
	}
	else if (function(a0) > function(a0 + h))
	{
		*(range) = a0;
		do
		{
			a0 += h;
		} while (function(a0)>function(a0 + h));
		*(range + 1) = a0 + h;
	}
	else
	{
		*(range + 1) = a0 + h;
		do
		{
			a0 -= h;
		} while (function(a0)<function(a0 + h));
		*(range) = a0 - h;
	}
}


/*************************************************
* @Description: 平分区间[a,b]为(n+1)等分
* @Input: a:区间左端点,b:区间右端点,A:等分点数组
* @Return: A[i]
*************************************************/
void findA(double a, double b, double *A) {
	for (int i = 0; i < N; i++)
	{
		*(A + i) = (b - a) / (N + 1)*(i + 1) + a;
	}
}

/*************************************************
* @Description: 求数组最小值序号
* @Input: A:数组
* @Return: 最小值序号
*************************************************/
int findMin(double *A) {
	double tem = *(A + 0);
	int num = 0;
	for (int i = 0; i < N; i++)
	{
		if (tem > *(A+i))
		{
			num = i;
			tem = *(A + i);
		}
	}
	return num;
}

/*************************************************
* @Description: 格点法求解主要函数
* @Input: a:区间左端点,b:区间右端点,varepsilon：精度
* @Return: A[i]
*************************************************/
double geDian(double a, double b, double varepsilon) {
	int num;
	double A[N], F[N];	//A为等分的点，F为等分点的函数值
	do
	{
		findA(a, b, A);	//得到等分点
		for (int i = 0; i < N; i++)	//得到等分点对应函数值
		{
			F[i] = function(A[i]);
		}
		num = findMin(F);
		switch (num)	//迭代
		{
		case 0:b = A[num + 1]; break;
		case N:a = A[num - 1]; break;
        case N - 1:a = A[(num - 1)]; break;
		default:a = A[num - 1]; b = A[num + 1]; break;
		}
	} while (b - a > 0 ? ((b - a) > varepsilon) : ((a - b) > varepsilon));
	return 1.0 / 2 * (a + b);
}

void main() {
	double a0 = A0, h = H, range[2] = { 0,0 }, varepsilon = V, result;
	jintui(a0, h, range);
	result = geDian(range[0],range[1],varepsilon);
	printf("Result:%f\n", result);
	system("pause");
}

2.3 过程分析

遇到的问题

2.3.1 C语言与C++的区别

我的循环是这样写的：

1
2
3

for (int i = 0; i < length; i++){
    //Do Some Thing
}

C++及Java等语言是支持这个的

但是我发现在【VC6.0】中的标准C语言不支持临时变量在for循环中定义

解决方法

在函数的开始申明变量，之后再使用变量

int i;
for (i = 0; i < length; i++){
    //Do Some Thing
}

2.3.2 数组的越界

当输入为一个很大的负数时，输出结果就会出现一个比较奇怪的数据

解决方法

调试程序，分析异常数据，最后发现是在 geDian 函数中，switch 语句少考虑了一种情况，导致数组越界。

将 geDian 函数中：

switch (num)	//迭代
		{
		case 0:a = a; b = A[num + 1]; break;
		case N:b = b; a = A[num - 1]; break;
		default:a = A[num - 1]; b = A[num + 1]; break;
		}

改为：

switch (num)	//迭代
		{
		case 0:b = A[num + 1]; break;
		case N:a = A[num - 1]; break;
        case N - 1:a = A[(num - 1)]; break;
		default:a = A[num - 1]; b = A[num + 1]; break;
		}

成功解决问题

2.4 结果分析

给定函数：

$f(\alpha)=\alpha^2+2\alpha+1$

其理论极小点为 $\alpha^*=-1$

给定初值 $\alpha_0=999$ ，给定步长 $h=10$

结果：

所求极小点 $\alpha^*=-1.000000$ ，满足精度，通过验证

3. 求解多元函数的最优解

实验日期：2018 年 11 月 22 日，2018 年 11 月 29 日

3.1 实验目的

加深对多元函数无约束优化方法的基本理论和算法步骤的理解。
培养学生独立编制计算机程序的能力。
掌握常用多元函数无约束优化方法程序的使用方法。

3.2 上机内容

3.2.1 算法概述

此处使用改进鲍威尔法进行求解。改进鲍威尔法用于求解多元函数函数极小点 $\alpha^*$ 。

其算法为：在每一轮迭代中总有个始点（第一轮的始点是任选的初始点）和 $n$ 个线性独立的搜索方向。从始点出发顺次沿 $n$ 个方向作一维搜索得一终点，由始点和终点决定了一个新的搜索方向。判断是否需要替换，若需要替换，用这个方向替换原来 $n$ 个方向中的一个，形成新的搜索方向组。这样就形成算法的循环。

替换的原则是去掉原方向组的第一个方向而将新方向排在原方向的最后。此外规定，从这一轮的搜索终点出发沿新的搜索方向作一维搜索而得到的极小点，作为下一轮迭代的始点。

3.2.2 程序框图

改进鲍威尔法

3.2.3 源程序

/*************************************************
* @Author: Fu_Qingchen
* @Date: 2018-11-22
* @Description: 改进鲍威尔法无约束优化
*************************************************/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>

#define A0	0	//给定初值1
#define A1	0	//给定初值2
#define V	0.000001	//给定精度

#define N	2	//函数变量个数
#define H	10	//格点法步长
#define N_1d	10	//格点法等分数-1
#define function_main(x1,x2) ((x1)*(x1)+2*(x2)*(x2)-4*(x1)-2*(x1)*(x2))	//要求解的函数

/*************************************************
* @Description: 进退法主函数
* @Input: a0:初始值,h:步长,*range:区间左右端点,*x:无约束问题的点，*S:搜索方向
* @Return: 区间左右端点
*************************************************/
void jintui(double a0, double h, double *range,double *x,double *S) {
	if (function_main(*(x)+a0*S[0],*(x+1)+a0*S[1]) == function_main(*(x)+(a0+h) * S[0], *(x + 1) + (a0+h) * S[1]))
	{
		*(range) = a0;
		*(range + 1) = a0 + h;
	}
	else if (function_main(*(x)+a0 * S[0], *(x + 1) + a0 * S[1]) > function_main(*(x)+(a0 + h) * S[0], *(x + 1) + (a0 + h) * S[1]))
	{
		*(range) = a0;
		do
		{
			a0 += h;
		} while (function_main(*(x)+a0 * S[0], *(x + 1) + a0 * S[1])>function_main(*(x)+(a0 + h) * S[0], *(x + 1) + (a0 + h) * S[1]));
		*(range + 1) = a0 + h;
	}
	else
	{
		*(range + 1) = a0 + h;
		do
		{
			a0 -= h;
		} while (function_main(*(x)+a0 * S[0], *(x + 1) + a0 * S[1])<function_main(*(x)+(a0 + h) * S[0], *(x + 1) + (a0 + h) * S[1]));
		*(range) = a0 - h;
	}
}

/*************************************************
* @Description: 平分区间[a,b]为(n+1)等分
* @Input: a:区间左端点,b:区间右端点,A:等分点数组
* @Return: A[i]
*************************************************/
void findA(double a, double b, double *A) {
	for (int i = 0; i < N_1d; i++)
	{
		*(A + i) = (b - a) / (N_1d + 1)*(i + 1) + a;
	}
}

/*************************************************
* @Description: 求数组最小值序号
* @Input: A:数组
* @Return: 最小值序号
*************************************************/
int findMin(double *A) {
	double tem = *(A + 0);
	int num = 0;
	for (int i = 0; i < N_1d; i++)
	{
		if (tem > *(A + i))
		{
			num = i;
			tem = *(A + i);
		}
	}
	return num;
}

/*************************************************
* @Description: 求数组最大值序号
* @Input: A:数组
* @Return: 最小值序号
*************************************************/
int findMax(double *A) {
	double tem = *(A + 0);
	int num = 0;
	for (int i = 0; i < N; i++)
	{
		if (tem < *(A + i))
		{
			num = i;
			tem = *(A + i);
		}
	}
	return num;
}

/*************************************************
* @Description: 格点法求解主要函数
* @Input: a:区间左端点,b:区间右端点,varepsilon:精度,*x:无约束问题的点，*S:搜索方向
* @Return: A[i]
*************************************************/
double geDian(double a, double b, double varepsilon,double *x, double *S) {
	int num;
	double A[N_1d], F[N_1d];	//A为等分的点，F为等分点的函数值
	do
	{
		findA(a, b, A);	//得到等分点
		for (int i = 0; i < N_1d; i++)	//得到等分点对应函数值
		{
			F[i] = function_main(*(x)+A[i] * S[0], *(x + 1) + A[i] * S[1]);
		}
		num = findMin(F);
		switch (num)	//迭代
		{
		case 0:b = A[(num + 1)]; break;
		case N_1d:a = A[(num - 1)]; break;
		case N_1d - 1:a = A[(num - 1)]; break;
		default:a = A[(num - 1)]; b = A[(num + 1)]; break;
		}
	} while (b - a > 0 ? ((b - a) > varepsilon) : ((a - b) > varepsilon));
	return 1.0 / 2 * (a + b);
}

/*************************************************
* @Description: 鲍威尔法求解主要函数
* @Input: *a0:初始点,varepsilon：精度
* @Return: *a0极小点
*************************************************/
int powell(double *a0, double varepsilon) {
	//设置初值
	double S[N+1][N+1] = { 0 };	//前进方向
	double S_new[N+1] = { 0 },sn1 = 0,sn2 = 0;	//构造的新方向
	double f[N + 1] = { 0 };	//迭代中的极小点
	double size[N] = { 0 };	//步长
	double delta[N] = { 0 };	//每前进一步的下降量
	int delta_max = 0;	//最大下降量对应序号
	double a00[N] = { 0 };	//初始点
	double a000[N] = { 0 };	//反射点
	double f1 = 0, f2 = 0, f3 = 0;	//方向有效性中的f1,f2,f3
	double h = H, range[2] = { 0,0 };	//一维搜索相关参数
	int i = 0, j = 0, k = 0;
	for (i = 0; i < N; i++)
	{
		S[i][i] = 1;
	}
	do
	{
		k++;
		printf("%d\n", k);
		f[0] = function_main(*(a0), *(a0 + 1));
		a00[0] = *(a0);
		a00[1] = *(a0 + 1);
		for (i = 0; i < N; i++)	//一维搜索
		{
			jintui(a0[i], h, range, a0, S[i]);
			size[i] = geDian(range[0], range[1], varepsilon, a0, S[i]);	//获取最优步长
			*(a0 + 0) = *(a0 + 0) + size[i] * S[i][0];	//获取极值点
			*(a0 + 1) = *(a0 + 1) + size[i] * S[i][1];
			f[i + 1] = function_main(*(a0), *(a0 + 1));	//获取极值点函数值
			delta[i] = f[i] - f[i + 1];
		}
		if (sqrt(pow(*(a0)-a00[0], 2) + pow(*(a0+1)-a00[1], 2))<varepsilon)	//收敛判断
		{
			return 0;
		}
		delta_max = findMax(delta);	//找到最大下降值
		for (i = 0; i < N; i++)
		{
			S_new[i] = *(a0 + i) - a00[i];	//构造新方向
			a000[i] = 2 * (*(a0 + i)) - a00[i];
		}
		sn1 = S_new[0];
		sn2 = S_new[1];
		f1 = f[0], f2 = f[N], f3 = function_main(a000[0], a000[1]);
		if (f3<f1 && ((f1 - 2 * f2 + f3)*pow(f1 - f2 - delta[delta_max], 2)<0.5*delta[delta_max] * pow(f1 - f3, 2)))	//方向有效
		{
			jintui(a0[0], h, range, a00, S_new);
			size[0] = geDian(range[0], range[1], varepsilon, a00, S_new);	//获取最优步长
			*(a0 + 0) = *(a00 + 0) + size[0] * S_new[0];	//获取极值点
			*(a0 + 1) = *(a00 + 1) + size[0] * S_new[1];
			f[0] = function_main(*(a0), *(a0 + 1));	//获取极值点函数值
			for (i = 0; i < N; i++)
			{
				if (i>=delta_max)
				{
					for ( j = 0; j < N; j++)
					{
						S[i - 1][j] = S[i][j];
					}
				}
			}
			S[N - 1][0] = sn1;
			S[N - 1][1] = sn2;
		}
		else	//方向无效
		{
			if (f3<f2)
			{
				*(a0 + 0) = a000[0];	//获取新初值
				*(a0 + 1) = a000[1];
			}
		}
	} while (1);
}

void main() {
	double a0[N], varepsilon = V;
	a0[0] = A0;
	a0[1] = A1;
	powell(a0, varepsilon);
	printf("x1:%lf\nx2:%lf\n", a0[0], a0[1]);
	system("pause");
}

3.3 过程分析

遇到的问题

3.3.1 参数的传递

在改进鲍威尔法中进行方向替换时，发现有一处一直替换不了

调试程序：powell 函数

printf("%lf;%lf\n",S_new[0],S_new[1]);
for (i = 0; i < N; i++)
{
	if (i>=delta_max)
	{
		for ( j = 0; j < N; j++)
		{
			S[i - 1][j] = S[i][j];
		}
	}
}
printf("%lf;%lf\n",S_new[0],S_new[1]);
S[N - 1][0] = S_new[0];
S[N - 1][1] = S_new[1];

发现两次输出不一样

解决方法

引入变量 sn1 和 sn2 通过它们来传递新方向

double sn1 = 0,sn2 = 0;	//引入的新方向
...
sn1 = S_new[0];
sn2 = S_new[1];
...
for (i = 0; i < N; i++)
{
	if (i>=delta_max)
	{
		for ( j = 0; j < N; j++)
		{
			S[i - 1][j] = S[i][j];
		}
	}
}
S[N - 1][0] = sn1;
S[N - 1][1] = sn2;