光滑牛顿法信赖域子问题求解python程序

以下是一个简单的Python程序，用于解决信赖域子问题：

import numpy as np
from scipy.optimize import minimize_scalar

def trust_region_subproblem(grad_f, hess_f, radius):
    """
    Trust region subproblem solver using the dogleg method
    """
    g = grad_f.reshape(-1, 1)
    B = hess_f
    p_u = -np.linalg.solve(B, g)
    norm_p_u = np.linalg.norm(p_u)

    if norm_p_u <= radius:
        # Return the unconstrained solution
        return p_u.flatten()

    p_b = -(g.T @ g) / (g.T @ B @ g) * g.flatten()
    norm_p_b = np.linalg.norm(p_b)

    if norm_p_b >= radius:
        # Return the boundary solution
        return radius / norm_p_u * p_u.flatten()

    # Solve the quadratic model
    a = np.linalg.norm(p_u - p_b) ** 2
    b = 2 * p_b.T @ (p_u - p_b)
    c = np.linalg.norm(p_b) ** 2 - radius ** 2
    alpha = (-b + np.sqrt(b ** 2 - 4 * a * c)) / (2 * a)
    p_d = p_b + alpha * (p_u - p_b)

    return p_d.flatten()

# Example usage
def f(x):
    return x[0] ** 2 + x[1] ** 2

def grad_f(x):
    return np.array([2 * x[0], 2 * x[1]])

def hess_f(x):
    return np.array([[2, 0], [0, 2]])

x0 = np.array([1, 1])
radius = 1.0

p = trust_region_subproblem(grad_f(x0), hess_f(x0), radius)
print(p)

在此示例中，我们定义了一个简单的二次函数f（x）= x [0] ** 2 + x [1] ** 2，并使用牛顿法求解其最小值。我们使用了Scipy库中的minimize_scalar函数来解决一维子问题。我们还定义了一个trust_region_subproblem函数，该函数使用dogleg方法解决了信赖域子问题。