Nonlinearly Constrained and Service Routine Examples¶
Example Code¶
NNLPF — Nonlinearly constrained minimization:
"""IMSL NNLPF example: nonlinearly constrained minimization.
Reproduces the IMSL NNLPF published example:
- Objective: f(x) = (x1-1)^2 + (x2-2)^2 + (x3-3)^2
- Nonlinear inequality constraint: x1^2 + x2^2 + x3^2 <= 14
(expressed as g(x) = 14 - x1^2 - x2^2 - x3^2 >= 0)
- Bounds: 0 <= xi <= 5 for i = 1, 2, 3
- Starting point: x0 = (0.5, 0.5, 0.5)
- Expected solution: x ≈ (1, 2, 3), f ≈ 0
Outputs:
- Table printed to stdout
- SVG plot saved to test_output/demo_imsl_nnlpf.svg
"""
from __future__ import annotations
from pathlib import Path
from typing import Dict
import matplotlib.pyplot as plt
import numpy as np
from optimization import minimize_nonlinear_constraints_fd, NonlinearConstraint
def run_demo_imsl_nnlpf() -> Dict[str, object]:
"""Run IMSL NNLPF example: nonlinearly constrained minimization.
Minimizes a sum-of-squares objective subject to a spherical inequality
constraint and box bounds, using finite-difference SLSQP (IMSL NNLPF).
Args:
None
Returns:
Dict[str, object]: Result dict with keys ``x``, ``fval``, ``n_iter``,
``n_fev``, and ``plot_path``.
"""
def f(x: np.ndarray) -> float:
"""Objective: (x1-1)^2 + (x2-2)^2 + (x3-3)^2.
Args:
x (np.ndarray): Decision variables of shape (3,).
Returns:
float: Objective value.
"""
return (x[0] - 1.0) ** 2 + (x[1] - 2.0) ** 2 + (x[2] - 3.0) ** 2
def g(x: np.ndarray) -> float:
"""Inequality constraint: 14 - x1^2 - x2^2 - x3^2 >= 0.
Args:
x (np.ndarray): Decision variables of shape (3,).
Returns:
float: Constraint residual (>= 0 means feasible).
"""
return 14.0 - x[0] ** 2 - x[1] ** 2 - x[2] ** 2
x0 = np.array([0.5, 0.5, 0.5])
x_expected = np.array([1.0, 2.0, 3.0])
f_expected = 0.0
result = minimize_nonlinear_constraints_fd(
f,
m=1,
me=0,
ibtype=0,
xlb=np.array([0.0, 0.0, 0.0]),
xub=np.array([5.0, 5.0, 5.0]),
constraints=[NonlinearConstraint(fun=g, type="ineq")],
x_guess=x0,
max_iter=500,
)
print("\nIMSL NNLPF Example: Nonlinearly Constrained Minimization")
print("-" * 60)
print(f"{'Component':<15} {'Solution':>15} {'Expected':>15}")
print("-" * 60)
for i, (xi, xe) in enumerate(zip(result.x, x_expected)):
print(f"{'x' + str(i + 1):<15} {xi:>15.6f} {xe:>15.6f}")
print("-" * 60)
print(f"{'f(x)':<15} {result.fval:>15.6f} {f_expected:>15.6f}")
print(f"{'Iterations':<15} {result.n_iter:>15}")
print(f"{'Func evals':<15} {result.n_fev:>15}")
print(f"{'Converged':<15} {str(result.success):>15}")
print("-" * 60)
output_dir = Path("test_output")
output_dir.mkdir(parents=True, exist_ok=True)
plot_path = output_dir / "demo_imsl_nnlpf.svg"
components = ["x1", "x2", "x3"]
x_sol = result.x
bar_w = 0.35
idx = np.arange(3)
fig, ax = plt.subplots(figsize=(7, 5))
ax.bar(idx, x_expected, bar_w, label="Expected (1, 2, 3)", color="#1f77b4", alpha=0.7)
ax.bar(idx + bar_w, x_sol, bar_w, label="NNLPF Solution", color="#d62728", alpha=0.7)
ax.set_xticks(idx + bar_w / 2)
ax.set_xticklabels(components)
ax.set_ylabel("Value")
ax.set_title(f"IMSL NNLPF: Nonlinear Constrained Solution (f={result.fval:.6f})")
ax.legend()
ax.grid(True, alpha=0.3, axis="y")
fig.tight_layout()
fig.savefig(plot_path, format="svg")
plt.close(fig)
return {
"x": result.x,
"fval": result.fval,
"n_iter": result.n_iter,
"n_fev": result.n_fev,
"plot_path": str(plot_path),
}
if __name__ == "__main__":
run_demo_imsl_nnlpf()
Service routines — CDGRD/FDGRD/CHGRD etc.:
"""IMSL service routines example: finite difference utilities.
Demonstrates all service routines for the IMSL Optimization chapter:
- central_difference_gradient (CDGRD)
- forward_difference_gradient (FDGRD)
- forward_difference_hessian (FDHES)
- gradient_difference_hessian (GDHES)
- check_gradient (CHGRD)
- check_hessian (CHHES)
- generate_starting_points (GGUES)
Uses Rosenbrock function: f(x) = 100*(x2-x1^2)^2 + (1-x1)^2
Outputs:
- Table printed to stdout
- SVG plot saved to test_output/demo_imsl_service.svg
"""
from __future__ import annotations
from pathlib import Path
from typing import Dict
import matplotlib.pyplot as plt
import numpy as np
from optimization import (
central_difference_gradient,
check_gradient,
check_hessian,
forward_difference_gradient,
forward_difference_hessian,
generate_starting_points,
gradient_difference_hessian,
)
def _rosenbrock(x: np.ndarray) -> float:
"""Rosenbrock function.
Args:
x (np.ndarray): Input vector of shape (2,).
Returns:
float: f(x) = 100*(x2-x1^2)^2 + (1-x1)^2.
"""
return 100.0 * (x[1] - x[0] ** 2) ** 2 + (1.0 - x[0]) ** 2
def _grad_rosenbrock(x: np.ndarray) -> np.ndarray:
"""Analytic gradient of Rosenbrock function.
Args:
x (np.ndarray): Input vector of shape (2,).
Returns:
np.ndarray: Gradient of shape (2,).
"""
return np.array([
-400.0 * (x[1] - x[0] ** 2) * x[0] - 2.0 * (1.0 - x[0]),
200.0 * (x[1] - x[0] ** 2),
])
def _hess_rosenbrock(x: np.ndarray) -> np.ndarray:
"""Analytic Hessian of Rosenbrock function.
Args:
x (np.ndarray): Input vector of shape (2,).
Returns:
np.ndarray: Hessian of shape (2, 2).
"""
return np.array([
[-400.0 * (x[1] - 3.0 * x[0] ** 2) + 2.0, -400.0 * x[0]],
[-400.0 * x[0], 200.0],
])
def run_demo_imsl_service() -> Dict[str, object]:
"""Run IMSL service routines example.
Args:
None
Returns:
Dict[str, object]: Result dict with keys ``cd_grad``, ``fd_grad``,
``fd_hess``, ``grad_check``, ``hess_check``, ``starting_points``,
and ``plot_path``.
"""
x0 = np.array([-1.2, 1.0])
cd_grad = central_difference_gradient(_rosenbrock, x0)
fd_grad = forward_difference_gradient(_rosenbrock, x0)
an_grad = _grad_rosenbrock(x0)
fd_hess = forward_difference_hessian(_rosenbrock, x0)
an_hess = _hess_rosenbrock(x0)
gd_hess = gradient_difference_hessian(_rosenbrock, _grad_rosenbrock, x0)
grad_check = check_gradient(_rosenbrock, _grad_rosenbrock, x0, tol=1e-3)
hess_check = check_hessian(_rosenbrock, _hess_rosenbrock, x0, tol=1e-3)
pts = generate_starting_points(n=2, n_points=10, bounds=[(-5.0, 5.0), (-5.0, 5.0)], seed=42)
print("\nIMSL Service Routines Example (Rosenbrock at x=(-1.2, 1.0))")
print("=" * 65)
print("\n1. Gradient approximations:")
print(f" Analytic: {an_grad}")
print(f" Central difference: {cd_grad}")
print(f" Forward difference: {fd_grad}")
print("\n2. Hessian approximations:")
print(" Analytic:")
print(f" {an_hess[0]}")
print(f" {an_hess[1]}")
print(" Forward difference:")
print(f" {fd_hess[0]}")
print(f" {fd_hess[1]}")
print(" Gradient-based:")
print(f" {gd_hess[0]}")
print(f" {gd_hess[1]}")
print("\n3. Gradient check (CHGRD):")
print(f" max_error = {grad_check['max_error']:.3e}, passed = {grad_check['passed']}")
print("\n4. Hessian check (CHHES):")
print(f" max_error = {hess_check['max_error']:.3e}, passed = {hess_check['passed']}")
print("\n5. Generated starting points (10 pts in [-5,5]^2):")
print(f" Shape: {pts.shape}")
print(f" Min per dim: {pts.min(axis=0)}")
print(f" Max per dim: {pts.max(axis=0)}")
output_dir = Path("test_output")
output_dir.mkdir(parents=True, exist_ok=True)
plot_path = output_dir / "demo_imsl_service.svg"
x1v = np.linspace(-3, 3, 200)
x2v = np.linspace(-2, 4, 200)
X1, X2 = np.meshgrid(x1v, x2v)
Z = np.vectorize(lambda a, b: _rosenbrock(np.array([a, b])))(X1, X2)
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].contourf(X1, X2, np.log1p(Z), levels=25, cmap="viridis")
axes[0].scatter(pts[:, 0], pts[:, 1], color="red", s=40, zorder=5, label="Starting points")
axes[0].scatter(*x0, color="white", s=80, marker="^", zorder=5, label="Eval point")
axes[0].set_title("Starting Points on Rosenbrock Landscape")
axes[0].set_xlabel("x1")
axes[0].set_ylabel("x2")
axes[0].legend(fontsize=8)
grad_labels = ["∂f/∂x1", "∂f/∂x2"]
grad_data = [an_grad, cd_grad, fd_grad]
grad_names = ["Analytic", "Central diff", "Forward diff"]
x_pos = np.arange(2)
bar_w = 0.25
colors = ["#1f77b4", "#ff7f0e", "#2ca02c"]
for k, (gd, gn, col) in enumerate(zip(grad_data, grad_names, colors)):
axes[1].bar(x_pos + k * bar_w, gd, bar_w, label=gn, color=col, alpha=0.8)
axes[1].set_xticks(x_pos + bar_w)
axes[1].set_xticklabels(grad_labels)
axes[1].set_title("Gradient Approximations at x=(-1.2, 1.0)")
axes[1].set_ylabel("Gradient value")
axes[1].legend()
axes[1].grid(True, alpha=0.3, axis="y")
fig.tight_layout()
fig.savefig(plot_path, format="svg")
plt.close(fig)
return {
"cd_grad": cd_grad,
"fd_grad": fd_grad,
"fd_hess": fd_hess,
"grad_check": grad_check,
"hess_check": hess_check,
"starting_points": pts,
"plot_path": str(plot_path),
}
if __name__ == "__main__":
run_demo_imsl_service()
Input (Console)¶
Run the nonlinearly constrained and service scripts from the package root:
python examples/example_imsl_nnlpf.py
python examples/example_imsl_service.py
Plot Output¶
Generated SVG plots:
Output Console¶
NNLPF console output:
IMSL NNLPF Example: Nonlinearly Constrained Minimization
------------------------------------------------------------
Component Solution Expected
------------------------------------------------------------
x1 1.000000 1.000000
x2 2.000000 2.000000
x3 3.000000 3.000000
------------------------------------------------------------
f(x) 0.000000 0.000000
Iterations 3
Func evals 14
Converged True
------------------------------------------------------------
Service routines console output:
Hessian check failed: max_error=4.218e-01 > tol=1.000e-03
IMSL Service Routines Example (Rosenbrock at x=(-1.2, 1.0))
=================================================================
1. Gradient approximations:
Analytic: [-215.6 -88. ]
Central difference: [-215.60000002 -88. ]
Forward difference: [-215.5999893 -87.99999857]
2. Hessian approximations:
Analytic:
[1330. 480.]
[480. 200.]
Forward difference:
[1329.57815693 479.97070312]
[479.97070312 199.99999976]
Gradient-based:
[1329.99998093 479.9999996 ]
[479.9999996 200. ]
3. Gradient check (CHGRD):
max_error = 2.208e-08, passed = True
4. Hessian check (CHHES):
max_error = 4.218e-01, passed = False
5. Generated starting points (10 pts in [-5,5]^2):
Shape: (10, 2)
Min per dim: [-4.92897428 -4.30931841]
Max per dim: [4.84961882 4.9854121 ]