最优化大作业（二）： 常用无约束最优化方法

2023-12-12 14:08
文章标签 方法 作业 常用 无约束 最优化

本文主要是介绍最优化大作业（二）： 常用无约束最优化方法，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

  • 问题描述

对以下优化问题

                                                   \small minf\left ( x \right )=x_1^2+x_2^2+x_1x_2-10x_1-4x_2+60

选取初始点\small X_0=\left [ 0,0 \right ]^T,\varepsilon =10^{-2}，分别用以下方法求解

（1）最速下降法；

（2）Newton法或修正Newton法；

（3）共轭梯度法。

 

  • 基本原理

（1）最速下降法

 

图1  最速下降法流程图

（2）Newton法

图2  Newton法流程图

 

（3）共轭梯度法

 

图3  共轭梯度法流程图

 

  • 实验结果

（1）最速下降法

迭代

次数

1

2

3

4

5

6

7

8

9

梯度

模值

 

5.4210

1.6680

0.9532

0.2933

0.1676

0.0516

0.0295

0.0091

搜索

方向

 

[9.00,

3.00]

[-1.71,

5.14]

[1.58,

0.53]

[-0.30,

0.90]

[0.28,

0.09]

[-0.05,

0.16]

[0.05,

0.02]

[-0.01,

0.03]

当前

迭代点

(1.00,

1.00)

(7.43,

3.14)

(6.77,

5.12)

(7.90,

5.50)

(7.78,

5.85)

(7.98,

5.91)

(7.96,

5.97)

(8.00,

5.98)

(7.99,

6.00)

当前迭代点值

47.00

14.8571

9.2057

8.2120

8.0373

8.0066

8.0012

8.0002

8.0000

表1  最速下降法迭代过程

图4  最速下降法迭代过程图

 

 

（2）Newton法

迭代次数

1

2

梯度模值

 

0.0000

搜索方向

 

[7.00,5.00]

当前迭代点

(1.00,1.00)

(8.00,6.00)

当前迭代点值

47.00

8.0000

表2  Newton法迭代过程

图5  Newton法迭代过程图

 

（3）共轭梯度法

迭代次数

1

2

3

梯度模值

 

5.4210

0.0000

搜索方向

 

[9.00,3.00]

[1.22,6.12]

当前迭代点

(1.00,1.00)

(7.43,3.14)

(8.00,6.00)

当前迭代点值

47.00

14.8571

8.0000

表3  共轭梯度法迭代过程

 

 

图6  共轭梯度法迭代过程图

 

对比结果可得，三种算法均得到同一个极值点（8, 6）。

 

  • 代码展示
import matplotlib.pyplot as plt
from sympy import *
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
t = symbols('t')# 优化目标函数
def fun1():x1, x2 = symbols('x1 x2')y = np.power(x1, 2) + np.power(x2, 2) - x1*x2 -10 * x1 - 4 *x2 +60return ydef fun2(x1, x2):return np.power(x1, 2) + np.power(x2, 2) - x1 * x2 - 10 * x1 - 4 * x2 + 60# 计算当前梯度
def cal_gradient_cur(X_cur):x1, x2 = symbols('x1 x2')f = fun1()g = [diff(f, x1), diff(f, x2)]g[0] = g[0].evalf(subs={x1:X_cur[0], x2:X_cur[1]})g[1] = g[1].evalf(subs={x1:X_cur[0], x2:X_cur[1]})return np.array(g)# 计算lambda, X1: 上一轮迭代点, X2: 本次迭代点
def cal_lambda(X1, X2):g1 = np.array(cal_gradient_cur(X1))g2 = np.array(cal_gradient_cur(X2))g1_norm_2 = np.sum(g1**2, axis=0)g2_norm_2 = np.sum(g2**2, axis=0)lamda = g2_norm_2 / g1_norm_2return lamda# 更新迭代点X
def update_X(X, P):return np.array(X + t*P)# 更新迭代点X
def update_X_cur(X, t, P):return np.array(X + t*P)# 计算最优步长
def cal_best_t(X_cur):x1, x2 = symbols('x1 x2')f = fun1()f_t = f.subs({x1: X_cur[0], x2: X_cur[1]})return solve(diff(f_t, t), t)# 计算梯度模值
def cal_g_norm_cur(X):g_cur = np.array(cal_gradient_cur(X), dtype=np.float32)return np.sqrt(np.sum(g_cur**2, axis=0))def draw(X0):plt.figure()ax = plt.axes(projection='3d')xx = np.arange(-20, 20, 0.1)yy = np.arange(-20, 20, 0.1)x1, x2 = np.meshgrid(xx, yy)Z = fun2(x1, x2)ax.plot_surface(x1, x2, Z, cmap='rainbow', alpha=0.5)X = np.array([X0[0]])Y = np.array([X0[1]])X, Y = np.meshgrid(X, Y)Z = fun2(X, Y)print("初始点：(%0.2f,%0.2f,%0.2f)" % (X, Y, Z))ax.scatter(X, Y, Z, c='k', s=20)ax.set_xlabel('X')ax.set_ylabel('Y')ax.set_zlabel('Z')# ax.legend()# ax.contour(X,Y,Z, zdim='z',offset=-2，cmap='rainbow)   #等高线图，要设置offset，为Z的最小值return axclass C_gradient(object):def __init__(self, X0):self.X0 = X0# 更新搜索方向def cal_P(self, g_cur, P1, lamda):P = -1 * g_cur + lamda*P1return np.array(P)def search(self):X1 = self.X0g_norm_cur = cal_g_norm_cur(X1)  # 计算梯度模值count = 0result = []if(g_norm_cur <= 0.01):print("极值点为({:.2f},{:.2f})".format(X1[0], X1[1]))x1, x2 = symbols('x1 x2')f = fun1()min_value = f.evalf(subs={x1: X1[0], x2: X1[1]})print("极小值为{:.4f}".format(min_value))else:P = -1 * cal_gradient_cur(X1)  # 计算当前负梯度方向while True:X2 = update_X(X1, P)t_cur = cal_best_t(X2)X2 = update_X_cur(X1, t_cur, P)g_cur = cal_gradient_cur(X2)g_norm_cur = cal_g_norm_cur(X2)x1, x2 = symbols('x1 x2')f = fun1()min_value = f.evalf(subs={x1: X2[0], x2: X2[1]})result.append([float(X2[0]), float(X2[1]), float(min_value)])print("当前梯度模值为{:.4f}，搜索方向为[{:.2f},{:.2f}]".format(g_norm_cur, P[0], P[1]), end=" ")print("极值点为({:.2f},{:.2f}), 极小值为{:.4f}".format(result[count][0], result[count][1], result[count][2]))if(g_norm_cur <= 0.01):return np.array(result)else:lamda = cal_lambda(X1, X2)P = self.cal_P(g_cur, P, lamda)X1 = X2count += 1def C_gradient_main():print("当前搜索方法为共轭梯度法")X0 = np.array([1, 1])ax = draw(X0)cg = C_gradient(X0)cg.ax = axresult = cg.search()ax.scatter(np.array([result[:, 0]]), np.array([result[:, 1]]), np.array([result[:, 2]]), c='k', s=20)plt.show()class steepest_gradient(object):def __init__(self, X0):self.X0 = X0def search(self):X1 = self.X0result = []count = 0while True:P = -1 * cal_gradient_cur(X1)  # 计算当前负梯度方向X2 = update_X(X1, P)t_cur = cal_best_t(X2)X2 = update_X_cur(X1, t_cur, P)g_norm_cur = cal_g_norm_cur(X2)x1, x2 = symbols('x1 x2')f = fun1()min_value = f.evalf(subs={x1: X2[0], x2: X2[1]})result.append([float(X2[0]), float(X2[1]), float(min_value)])print("当前梯度模值为{:.4f}，搜索方向为[{:.2f},{:.2f}]".format(g_norm_cur, P[0], P[1]), end=" ")print("极值点为({:.2f},{:.2f}), 极小值为{:.4f}".format(result[count][0], result[count][1], result[count][2]))if(g_norm_cur <= 0.01):return np.array(result)else:X1 = X2count += 1def steepest_gradient_main():print("当前搜索方法为最速下降法")X0 = np.array([1, 1])ax = draw(X0)a = steepest_gradient(X0)a.ax = axresult = a.search()ax.scatter(np.array([result[:, 0]]), np.array([result[:, 1]]), np.array([result[:, 2]]), c='k', s=20)plt.show()class Newton(object):def __init__(self, X0):self.X0 = X0def cal_hesse(self):return np.array([[2, -1], [-1, 2]])def search(self):X1 = self.X0count = 0result = []while True:g = cal_gradient_cur(X1)g = g.reshape((1, 2)).Th = np.linalg.inv(self.cal_hesse())P = -1 * np.dot(h, g).ravel()X2 = update_X(X1, P)t_cur = cal_best_t(X2)X2 = update_X_cur(X1, t_cur, P)x1, x2 = symbols('x1 x2')f = fun1()min_value = f.evalf(subs={x1: X2[0], x2: X2[1]})g_norm_cur = cal_g_norm_cur(X2)result.append([float(X2[0]), float(X2[1]), float(min_value)])print("当前梯度模值为{:.4f}，搜索方向为[{:.2f},{:.2f}]".format(g_norm_cur, P[0], P[1]), end=" ")print("极值点为({:.2f},{:.2f}), 极小值为{:.4f}".format(result[count][0], result[count][1], result[count][2]))if(g_norm_cur <= 0.01):return np.array(result)else:X1 = X2count += 1def newton_main():print("当前搜索方法为newton法")X0 = np.array([1, 1])ax = draw(X0)b = Newton(X0)result = b.search()ax.scatter(np.array([result[:, 0]]), np.array([result[:, 1]]), np.array([result[:, 2]]), c='k', s=20)plt.show()if __name__ == '__main__':steepest_gradient_main()newton_main()C_gradient_main()

 

这篇关于最优化大作业（二）： 常用无约束最优化方法的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/484903。

相关文章

C# 比较两个list 之间元素差异的常用方法

C# 比较两个list 之间元素差异的常用方法

《C#比较两个list之间元素差异的常用方法》：本文主要介绍C#比较两个list之间元素差异,本文通过实例代码给大家介绍的非常详细,对大家的学习或工作具有一定的参考借鉴价值,需要的朋友参考下吧... 目录1. 使用Except方法2. 使用Except的逆操作3. 使用LINQ的Join,GroupJoin
阅读更多...
MySQL查询JSON数组字段包含特定字符串的方法

MySQL查询JSON数组字段包含特定字符串的方法

《MySQL查询JSON数组字段包含特定字符串的方法》在MySQL数据库中,当某个字段存储的是JSON数组,需要查询数组中包含特定字符串的记录时传统的LIKE语句无法直接使用,下面小编就为大家介绍两种... 目录问题背景解决方案对比1. 精确匹配方案（推荐）2. 模糊匹配方案参数化查询示例使用场景建议性能优
阅读更多...
关于集合与数组转换实现方法

关于集合与数组转换实现方法

《关于集合与数组转换实现方法》：本文主要介绍关于集合与数组转换实现方法,具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地方,望不吝赐教... 目录1、Arrays.asList()1.1、方法作用1.2、内部实现1.3、修改元素的影响1.4、注意事项2、list.toArray()2.1、方
阅读更多...
Python中注释使用方法举例详解

Python中注释使用方法举例详解

《Python中注释使用方法举例详解》在Python编程语言中注释是必不可少的一部分,它有助于提高代码的可读性和维护性,：本文主要介绍Python中注释使用方法的相关资料,需要的朋友可以参考下... 目录一、前言二、什么是注释？示例：三、单行注释语法：以 China编程# 开头，后面的内容为注释内容示例：示例：四
阅读更多...
一文详解Git中分支本地和远程删除的方法

一文详解Git中分支本地和远程删除的方法

《一文详解Git中分支本地和远程删除的方法》在使用Git进行版本控制的过程中,我们会创建多个分支来进行不同功能的开发,这就容易涉及到如何正确地删除本地分支和远程分支,下面我们就来看看相关的实现方法吧... 目录技术背景实现步骤删除本地分支删除远程www.chinasem.cn分支同步删除信息到其他机器示例步骤
阅读更多...
python常用的正则表达式及作用

python常用的正则表达式及作用

《python常用的正则表达式及作用》正则表达式是处理字符串的强大工具,Python通过re模块提供正则表达式支持,本文给大家介绍python常用的正则表达式及作用详解,感兴趣的朋友跟随小编一起看看吧... 目录python常用正则表达式及作用基本匹配模式常用正则表达式示例常用量词边界匹配分组和捕获常用re
阅读更多...
在Golang中实现定时任务的几种高效方法

在Golang中实现定时任务的几种高效方法

《在Golang中实现定时任务的几种高效方法》本文将详细介绍在Golang中实现定时任务的几种高效方法,包括time包中的Ticker和Timer、第三方库cron的使用,以及基于channel和go... 目录背景介绍目的和范围预期读者文档结构概述术语表核心概念与联系故事引入核心概念解释核心概念之间的关系
阅读更多...
在Linux终端中统计非二进制文件行数的实现方法

在Linux终端中统计非二进制文件行数的实现方法

《在Linux终端中统计非二进制文件行数的实现方法》在Linux系统中,有时需要统计非二进制文件（如CSV、TXT文件）的行数,而不希望手动打开文件进行查看,例如,在处理大型日志文件、数据文件时,了解... 目录在linux终端中统计非二进制文件的行数技术背景实现步骤1. 使用wc命令2. 使用grep命令
阅读更多...
Python中Tensorflow无法调用GPU问题的解决方法

Python中Tensorflow无法调用GPU问题的解决方法

《Python中Tensorflow无法调用GPU问题的解决方法》文章详解如何解决TensorFlow在Windows无法识别GPU的问题,需降级至2.10版本,安装匹配CUDA11.2和cuDNN... 当用以下代码查看GPU数量时，gpuspython返回的是一个空列表，说明tensorflow没有找到
阅读更多...
XML重复查询一条Sql语句的解决方法

XML重复查询一条Sql语句的解决方法

《XML重复查询一条Sql语句的解决方法》文章分析了XML重复查询与日志失效问题,指出因DTO缺少@Data注解导致日志无法格式化、空指针风险及参数穿透,进而引发性能灾难,解决方案为在Controll... 目录一、核心问题：从SQL重复执行到日志失效二、根因剖析：DTO断裂引发的级联故障三、解决方案：修复
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2