# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at https://mozilla.org/MPL/2.0/.
# SPDX-License-Identifier: MPL-2.0
import time
import numpy as np
from typing import Callable, Any
from VeraGridEngine.basic_structures import Vec
from VeraGridEngine.Utils.NumericalMethods.common import (ConvexMethodResult, ConvexFunctionResult,
check_function_and_args)
from VeraGridEngine.Utils.NumericalMethods.sparse_solve import get_linear_solver
from VeraGridEngine.basic_structures import Logger
linear_solver = get_linear_solver()
[docs]
def newton_raphson(func: Callable[[Vec, bool, Any], ConvexFunctionResult],
func_args: Any,
x0: Vec,
tol: float = 1e-6,
max_iter: int = 10,
trust: float = 1.0,
verbose: int = 0,
logger: Logger = Logger()) -> ConvexMethodResult:
"""
Newton-Raphson with Line search to solve:
min: error(g(x))
s.t.
g(x) = 0
:param func: function to optimize, it may or may not include the Jacobian matrix
the function must have x: Vec and calc_jacobian: bool as the first arguments
the function must return an instance of ConvexFunctionResult that contains
the function vector and optionally the derivative when calc_jacobian=True
:param func_args: Tuple of static arguments to call the evaluation function
:param x0: Array of initial solutions
:param tol: Error tolerance
:param max_iter: Maximum number of iterations
:param trust: trust amount in the derivative length correctness
:param verbose: Display console information
:param logger: Logger instance
:return: ConvexMethodResult
"""
start = time.time()
if not check_function_and_args(func, func_args, 2):
raise Exception(f'Invalid function arguments, required {", ".join(func.__code__.co_varnames)}')
# evaluation of the initial point
x = x0.copy()
ret = func(x, True, *func_args) # compute the Jacobian too
error = ret.compute_f_error()
converged = error < tol
iteration = 0
error_evolution = np.zeros(max_iter + 1)
trust0 = trust if trust <= 1.0 else 1.0 # trust radius in NR should not be greater than 1
# save the error evolution
error_evolution[iteration] = error
if verbose > 0:
print(f'It {iteration}, error {error}, converged {converged}, x {x}, dx not computed yet')
if converged:
return ConvexMethodResult(x=x0,
error=error,
converged=converged,
iterations=iteration,
elapsed=time.time() - start,
error_evolution=error_evolution)
else:
while not converged and iteration < max_iter:
print("(newton_raphson.py) Iteration: ", iteration)
# compute update step
try:
# compute update step: J x Ξx = Ξg
dx = linear_solver(ret.J, ret.f)
if np.isnan(dx).any():
logger.add_error(f"Newton-Raphson's Jacobian is singular @iter {iteration}:")
print("(newton_raphson.py) Singular Jacobian")
return ConvexMethodResult(x=x,
error=error,
converged=converged,
iterations=iteration,
elapsed=time.time() - start,
error_evolution=error_evolution)
except RuntimeError:
logger.add_error(f"Newton-Raphson's Jacobian is singular @iter {iteration}:")
print("(newton_raphson.py) Singular Jacobian")
return ConvexMethodResult(x=x,
error=error,
converged=converged,
iterations=iteration,
elapsed=time.time() - start,
error_evolution=error_evolution)
mu = trust0
back_track_condition = True
l_iter = 0
while back_track_condition and mu > tol:
x2 = x - mu * dx
ret2 = func(x2, False, *func_args) # do not compute the Jacobian
error2 = ret2.compute_f_error()
# change mu for the next iteration
mu *= 0.5 # acceleration_parameter
# keep back-tracking?
back_track_condition = error2 > error
l_iter += 1
if not back_track_condition:
# accept the solution
x = x2
if back_track_condition:
# this means that not even the backtracking was able to correct
# the solution, so terminate
logger.add_warning(f"Newton-Raphson's stagnated @iter {iteration}:")
return ConvexMethodResult(x=x,
error=error,
converged=converged,
iterations=iteration,
elapsed=time.time() - start,
error_evolution=error_evolution)
# update equalities increment and the jacobian
ret = func(x, True, *func_args)
# compute error
error = ret.compute_f_error()
# determine the convergence condition
converged = error <= tol
# update iteration counter
iteration += 1
# save the error evolution
error_evolution[iteration] = error
if verbose > 0:
print(f'It {iteration}, error {error}, converged {converged}, x {x}, dx {dx}')
return ConvexMethodResult(x=x,
error=error,
converged=converged,
iterations=iteration,
elapsed=time.time() - start,
error_evolution=error_evolution)