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

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

相关文章

Android 12解决push framework.jar无法开机的方法小结

《Android12解决pushframework.jar无法开机的方法小结》:本文主要介绍在Android12中解决pushframework.jar无法开机的方法,包括编译指令、框架层和s... 目录1. android 编译指令1.1 framework层的编译指令1.2 替换framework.ja

在.NET平台使用C#为PDF添加各种类型的表单域的方法

《在.NET平台使用C#为PDF添加各种类型的表单域的方法》在日常办公系统开发中,涉及PDF处理相关的开发时,生成可填写的PDF表单是一种常见需求,与静态PDF不同,带有**表单域的文档支持用户直接在... 目录引言使用 PdfTextBoxField 添加文本输入域使用 PdfComboBoxField

SQLyog中DELIMITER执行存储过程时出现前置缩进问题的解决方法

《SQLyog中DELIMITER执行存储过程时出现前置缩进问题的解决方法》在SQLyog中执行存储过程时出现的前置缩进问题,实际上反映了SQLyog对SQL语句解析的一个特殊行为,本文给大家介绍了详... 目录问题根源正确写法示例永久解决方案为什么命令行不受影响?最佳实践建议问题根源SQLyog的语句分

Python中pywin32 常用窗口操作的实现

《Python中pywin32常用窗口操作的实现》本文主要介绍了Python中pywin32常用窗口操作的实现,pywin32主要的作用是供Python开发者快速调用WindowsAPI的一个... 目录获取窗口句柄获取最前端窗口句柄获取指定坐标处的窗口根据窗口的完整标题匹配获取句柄根据窗口的类别匹配获取句

Java 中的 @SneakyThrows 注解使用方法(简化异常处理的利与弊)

《Java中的@SneakyThrows注解使用方法(简化异常处理的利与弊)》为了简化异常处理,Lombok提供了一个强大的注解@SneakyThrows,本文将详细介绍@SneakyThro... 目录1. @SneakyThrows 简介 1.1 什么是 Lombok?2. @SneakyThrows

判断PyTorch是GPU版还是CPU版的方法小结

《判断PyTorch是GPU版还是CPU版的方法小结》PyTorch作为当前最流行的深度学习框架之一,支持在CPU和GPU(NVIDIACUDA)上运行,所以对于深度学习开发者来说,正确识别PyTor... 目录前言为什么需要区分GPU和CPU版本?性能差异硬件要求如何检查PyTorch版本?方法1:使用命

Qt实现网络数据解析的方法总结

《Qt实现网络数据解析的方法总结》在Qt中解析网络数据通常涉及接收原始字节流,并将其转换为有意义的应用层数据,这篇文章为大家介绍了详细步骤和示例,感兴趣的小伙伴可以了解下... 目录1. 网络数据接收2. 缓冲区管理(处理粘包/拆包)3. 常见数据格式解析3.1 jsON解析3.2 XML解析3.3 自定义

SpringMVC 通过ajax 前后端数据交互的实现方法

《SpringMVC通过ajax前后端数据交互的实现方法》:本文主要介绍SpringMVC通过ajax前后端数据交互的实现方法,本文给大家介绍的非常详细,对大家的学习或工作具有一定的参考借鉴价... 在前端的开发过程中,经常在html页面通过AJAX进行前后端数据的交互,SpringMVC的controll

Java中的工具类命名方法

《Java中的工具类命名方法》:本文主要介绍Java中的工具类究竟如何命名,本文给大家介绍的非常详细,对大家的学习或工作具有一定的参考借鉴价值,需要的朋友参考下吧... 目录Java中的工具类究竟如何命名?先来几个例子几种命名方式的比较到底如何命名 ?总结Java中的工具类究竟如何命名?先来几个例子JD

Spring Security自定义身份认证的实现方法

《SpringSecurity自定义身份认证的实现方法》:本文主要介绍SpringSecurity自定义身份认证的实现方法,下面对SpringSecurity的这三种自定义身份认证进行详细讲解,... 目录1.内存身份认证(1)创建配置类(2)验证内存身份认证2.JDBC身份认证(1)数据准备 (2)配置依