最优化大作业(二): 常用无约束最优化方法

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

相关文章

Java学习手册之Filter和Listener使用方法

《Java学习手册之Filter和Listener使用方法》:本文主要介绍Java学习手册之Filter和Listener使用方法的相关资料,Filter是一种拦截器,可以在请求到达Servl... 目录一、Filter(过滤器)1. Filter 的工作原理2. Filter 的配置与使用二、Listen

Pandas统计每行数据中的空值的方法示例

《Pandas统计每行数据中的空值的方法示例》处理缺失数据(NaN值)是一个非常常见的问题,本文主要介绍了Pandas统计每行数据中的空值的方法示例,具有一定的参考价值,感兴趣的可以了解一下... 目录什么是空值?为什么要统计空值?准备工作创建示例数据统计每行空值数量进一步分析www.chinasem.cn处

Python的time模块一些常用功能(各种与时间相关的函数)

《Python的time模块一些常用功能(各种与时间相关的函数)》Python的time模块提供了各种与时间相关的函数,包括获取当前时间、处理时间间隔、执行时间测量等,:本文主要介绍Python的... 目录1. 获取当前时间2. 时间格式化3. 延时执行4. 时间戳运算5. 计算代码执行时间6. 转换为指

Windows 上如果忘记了 MySQL 密码 重置密码的两种方法

《Windows上如果忘记了MySQL密码重置密码的两种方法》:本文主要介绍Windows上如果忘记了MySQL密码重置密码的两种方法,本文通过两种方法结合实例代码给大家介绍的非常详细,感... 目录方法 1:以跳过权限验证模式启动 mysql 并重置密码方法 2:使用 my.ini 文件的临时配置在 Wi

MySQL重复数据处理的七种高效方法

《MySQL重复数据处理的七种高效方法》你是不是也曾遇到过这样的烦恼:明明系统测试时一切正常,上线后却频频出现重复数据,大批量导数据时,总有那么几条不听话的记录导致整个事务莫名回滚,今天,我就跟大家分... 目录1. 重复数据插入问题分析1.1 问题本质1.2 常见场景图2. 基础解决方案:使用异常捕获3.

最详细安装 PostgreSQL方法及常见问题解决

《最详细安装PostgreSQL方法及常见问题解决》:本文主要介绍最详细安装PostgreSQL方法及常见问题解决,介绍了在Windows系统上安装PostgreSQL及Linux系统上安装Po... 目录一、在 Windows 系统上安装 PostgreSQL1. 下载 PostgreSQL 安装包2.

Python正则表达式语法及re模块中的常用函数详解

《Python正则表达式语法及re模块中的常用函数详解》这篇文章主要给大家介绍了关于Python正则表达式语法及re模块中常用函数的相关资料,正则表达式是一种强大的字符串处理工具,可以用于匹配、切分、... 目录概念、作用和步骤语法re模块中的常用函数总结 概念、作用和步骤概念: 本身也是一个字符串,其中

SQL中redo log 刷⼊磁盘的常见方法

《SQL中redolog刷⼊磁盘的常见方法》本文主要介绍了SQL中redolog刷⼊磁盘的常见方法,将redolog刷入磁盘的方法确保了数据的持久性和一致性,下面就来具体介绍一下,感兴趣的可以了解... 目录Redo Log 刷入磁盘的方法Redo Log 刷入磁盘的过程代码示例(伪代码)在数据库系统中,r

usb接口驱动异常问题常用解决方案

《usb接口驱动异常问题常用解决方案》当遇到USB接口驱动异常时,可以通过多种方法来解决,其中主要就包括重装USB控制器、禁用USB选择性暂停设置、更新或安装新的主板驱动等... usb接口驱动异常怎么办,USB接口驱动异常是常见问题,通常由驱动损坏、系统更新冲突、硬件故障或电源管理设置导致。以下是常用解决

Python实现图片分割的多种方法总结

《Python实现图片分割的多种方法总结》图片分割是图像处理中的一个重要任务,它的目标是将图像划分为多个区域或者对象,本文为大家整理了一些常用的分割方法,大家可以根据需求自行选择... 目录1. 基于传统图像处理的分割方法(1) 使用固定阈值分割图片(2) 自适应阈值分割(3) 使用图像边缘检测分割(4)