逐次二次計画法（SQP Methods）

概要

逐次二次計画法（Sequential Quadratic Programming, SQP）は、制約付き非線形最適化問題に対する最も効果的な手法の一つです。各反復で二次計画問題を解くことにより、元の非線形問題を近似的に解き、ニュートン法と同様の二次収束性を実現します。

問題設定

一般的な制約付き非線形最適化問題：

$$\begin{align} \min_{x \in \mathbb{R}^n} &\quad f(x) \\ \text{subject to} &\quad h_j(x) = 0, \quad j = 1, 2, \ldots, p \\ &\quad g_i(x) \leq 0, \quad i = 1, 2, \ldots, m \end{align}$$

ここで：

• $f: \mathbb{R}^n \to \mathbb{R}$ は目的関数
• $h_j: \mathbb{R}^n \to \mathbb{R}$ は等式制約関数
• $g_i: \mathbb{R}^n \to \mathbb{R}$ は不等式制約関数

KKT条件と基本理論

ラグランジュ関数

$$L(x, \lambda, \mu) = f(x) + \sum_{j=1}^p \lambda_j h_j(x) + \sum_{i=1}^m \mu_i g_i(x)$$

KKT条件

最適解 $(x^*, \lambda^*, \mu^*)$ では以下が成立：

1. 停止性条件: $\nabla_x L(x^*, \lambda^*, \mu^*) = 0$
2. 実行可能性: $h_j(x^*) = 0$, $g_i(x^*) \leq 0$
3. 双対実行可能性: $\mu_i^* \geq 0$
4. 相補スラック性: $\mu_i^* g_i(x^*) = 0$

SQPの基本アイデア

ニュートン法との類推

無制約最適化のニュートン法：

$$x_{k+1} = x_k - [\nabla^2 f(x_k)]^{-1} \nabla f(x_k)$$

制約付き問題では、KKT条件にニュートン法を適用：

$$\begin{bmatrix} \nabla^2_x L_k & A_k^T \\ A_k & 0 \end{bmatrix} \begin{bmatrix} d_k \\ \lambda_{k+1} \end{bmatrix} = -\begin{bmatrix} \nabla_x L_k \\ h(x_k) \end{bmatrix}$$

ここで、$A_k = \nabla h(x_k)^T$ は制約のヤコビアン行列です。

SQP部分問題

各反復 $k$ で以下の二次計画問題を解きます：

$$\begin{align} \min_{d \in \mathbb{R}^n} &\quad \frac{1}{2} d^T B_k d + \nabla f(x_k)^T d \\ \text{subject to} &\quad \nabla h_j(x_k)^T d + h_j(x_k) = 0, \quad j = 1, \ldots, p \\ &\quad \nabla g_i(x_k)^T d + g_i(x_k) \leq 0, \quad i = 1, \ldots, m \end{align}$$

ここで：

• $d$ は探索方向
• $B_k$ はラグランジュ関数のヘッセ行列の近似（通常はBFGS更新）

アルゴリズム

基本SQPアルゴリズム

Algorithm: Basic SQP Method
1. 初期点 x₀, 初期ラグランジュ乗数 λ₀, μ₀ を選択
2. 初期ヘッセ行列近似 B₀ = I を設定
3. for k = 0, 1, 2, ... do:
   4. SQP部分問題を解いて (dₖ, λₖ₊₁, μₖ₊₁) を求める
   5. Merit関数に基づく直線探索でステップサイズ αₖ を決定
   6. 更新: xₖ₊₁ = xₖ + αₖdₖ
   7. ヘッセ行列近似を更新: Bₖ₊₁ = BFGS_update(Bₖ, sₖ, yₖ)
   8. 収束判定

Merit関数

SQPでは、目的関数と制約違反を組み合わせたメリット関数を使用：

L1メリット関数

$$\phi_1(x, \sigma) = f(x) + \sigma \left( \sum_{j=1}^p |h_j(x)| + \sum_{i=1}^m \max(0, g_i(x)) \right)$$

L2メリット関数

$$\phi_2(x, \sigma) = f(x) + \sigma \left( \sum_{j=1}^p h_j(x)^2 + \sum_{i=1}^m [\max(0, g_i(x))]^2 \right)$$

拡張ラグランジュメリット関数

$$\phi_{AL}(x, \lambda, \mu, \sigma) = f(x) + \sum_{j=1}^p \lambda_j h_j(x) + \sum_{i=1}^m \mu_i \max(0, g_i(x)) + \frac{\sigma}{2} \left( \sum_{j=1}^p h_j(x)^2 + \sum_{i=1}^m [\max(0, g_i(x))]^2 \right)$$

実装例

Python実装

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from scipy.sparse import csc_matrix
import warnings

class SQPOptimizer:
    """SQP法の実装クラス"""
    
    def __init__(self, merit_function='l1', hessian_update='bfgs'):
        """
        Parameters:
        merit_function: 'l1', 'l2', 'augmented_lagrangian'
        hessian_update: 'bfgs', 'exact'
        """
        self.merit_function = merit_function
        self.hessian_update = hessian_update
    
    def optimize(self, f, grad_f, h=None, grad_h=None, g=None, grad_g=None,
                x0=None, lambda0=None, mu0=None, max_iter=100, tol=1e-6, 
                sigma=1.0, verbose=False):
        """
        SQP法による制約付き最適化
        
        Parameters:
        f: 目的関数
        grad_f: 目的関数の勾配
        h: 等式制約関数のリスト
        grad_h: 等式制約の勾配のリスト
        g: 不等式制約関数のリスト  
        grad_g: 不等式制約の勾配のリスト
        """
        
        # デフォルト値の設定
        if h is None:
            h, grad_h = [], []
        if g is None:
            g, grad_g = [], []
        if x0 is None:
            x0 = np.zeros(2)
        if lambda0 is None:
            lambda0 = np.zeros(len(h))
        if mu0 is None:
            mu0 = np.zeros(len(g))
        
        n = len(x0)
        p = len(h)
        m = len(g)
        
        # 初期化
        x = x0.copy()
        lambda_eq = lambda0.copy()
        mu_ineq = mu0.copy()
        B = np.eye(n)  # 初期ヘッセ行列近似
        
        history = {
            'x_values': [x.copy()],
            'f_values': [f(x)],
            'constraint_violations': [],
            'kkt_residuals': [],
            'merit_values': [],
            'step_sizes': []
        }
        
        if verbose:
            print(f"SQP最適化開始 ({self.merit_function} メリット関数)")
            print(f"初期点: {x}")
            print(f"初期目的関数値: {f(x):.6f}")
        
        for k in range(max_iter):
            # 現在の制約値と勾配
            h_vals = np.array([hj(x) for hj in h]) if h else np.array([])
            g_vals = np.array([gi(x) for gi in g]) if g else np.array([])
            
            # 制約のヤコビアン
            if h:
                A_eq = np.array([grad_hj(x) for grad_hj in grad_h])
            else:
                A_eq = np.zeros((0, n))
            
            if g:
                A_ineq = np.array([grad_gi(x) for grad_gi in grad_g])
            else:
                A_ineq = np.zeros((0, n))
            
            # 制約違反の計算
            eq_violation = np.linalg.norm(h_vals) if len(h_vals) > 0 else 0
            ineq_violation = np.sum(np.maximum(0, g_vals)) if len(g_vals) > 0 else 0
            total_violation = eq_violation + ineq_violation
            
            # KKT残差の計算
            grad_lag = self._compute_lagrangian_gradient(
                x, grad_f, A_eq, A_ineq, lambda_eq, mu_ineq)
            kkt_residual = np.linalg.norm(grad_lag)
            
            # 履歴の更新
            history['constraint_violations'].append(total_violation)
            history['kkt_residuals'].append(kkt_residual)
            
            # メリット関数値の計算
            merit_val = self._compute_merit_function(
                x, f, h_vals, g_vals, lambda_eq, mu_ineq, sigma)
            history['merit_values'].append(merit_val)
            
            if verbose and k % 5 == 0:
                print(f"反復 {k}: f = {f(x):.6f}, 制約違反 = {total_violation:.6f}, "
                      f"KKT残差 = {kkt_residual:.6f}")
            
            # 収束判定
            if kkt_residual < tol and total_violation < tol:
                if verbose:
                    print(f"\n収束しました (反復回数: {k})")
                break
            
            # SQP部分問題を解く
            try:
                d, lambda_new, mu_new = self._solve_qp_subproblem(
                    x, grad_f(x), B, h, grad_h, g, grad_g)
                
                # Merit関数に基づく直線探索
                alpha = self._line_search_merit(
                    x, d, f, h, g, lambda_eq, mu_ineq, sigma)
                
                # 更新
                x_new = x + alpha * d
                
                # ラグランジュ乗数の更新
                if len(lambda_new) > 0:
                    lambda_eq = lambda_new
                if len(mu_new) > 0:
                    mu_ineq = np.maximum(0, mu_new)  # 非負性の確保
                
                # ヘッセ行列近似の更新
                if alpha > 1e-8:  # 十分なステップが取れた場合
                    s = x_new - x
                    y = self._compute_lagrangian_gradient_difference(
                        x, x_new, grad_f, grad_h, grad_g, lambda_eq, mu_ineq)
                    
                    if np.dot(s, y) > 1e-12:  # 曲率条件
                        B = self._update_hessian_approximation(B, s, y)
                
                x = x_new
                
                # 履歴の更新
                history['x_values'].append(x.copy())
                history['f_values'].append(f(x))
                history['step_sizes'].append(alpha)
                
            except Exception as e:
                if verbose:
                    print(f"反復 {k} でエラーが発生: {e}")
                break
        
        return x, lambda_eq, mu_ineq, history
    
    def _solve_qp_subproblem(self, x, grad_f_x, B, h, grad_h, g, grad_g):
        """SQP部分問題（二次計画問題）を解く"""
        
        n = len(x)
        
        # 制約の評価
        h_vals = np.array([hj(x) for hj in h]) if h else np.array([])
        g_vals = np.array([gi(x) for gi in g]) if g else np.array([])
        
        # 制約のヤコビアン
        A_eq = np.array([grad_hj(x) for grad_hj in grad_h]) if h else None
        A_ineq = np.array([grad_gi(x) for grad_gi in grad_g]) if g else None
        
        # scipyのminimizeを使用してQPを解く
        def qp_objective(d):
            return 0.5 * d.T @ B @ d + grad_f_x.T @ d
        
        def qp_objective_grad(d):
            return B @ d + grad_f_x
        
        # 制約の設定
        constraints = []
        
        # 等式制約
        if h:
            for i, (A_row, h_val) in enumerate(zip(A_eq, h_vals)):
                constraints.append({
                    'type': 'eq',
                    'fun': lambda d, A=A_row, b=h_val: A @ d + b,
                    'jac': lambda d, A=A_row: A
                })
        
        # 不等式制約
        if g:
            for i, (A_row, g_val) in enumerate(zip(A_ineq, g_vals)):
                constraints.append({
                    'type': 'ineq', 
                    'fun': lambda d, A=A_row, b=g_val: -(A @ d + b),
                    'jac': lambda d, A=A_row: -A
                })
        
        # QP問題を解く
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            result = minimize(qp_objective, np.zeros(n), method='SLSQP',
                            jac=qp_objective_grad, constraints=constraints,
                            options={'ftol': 1e-12, 'disp': False})
        
        d = result.x
        
        # ラグランジュ乗数の抽出（近似）
        if hasattr(result, 'v') and result.v is not None:
            # SLSQPの場合、result.vにラグランジュ乗数が格納される
            multipliers = result.v
            lambda_eq = multipliers[:len(h)] if h else np.array([])
            mu_ineq = multipliers[len(h):] if g else np.array([])
        else:
            # ラグランジュ乗数が取得できない場合は近似
            lambda_eq = np.zeros(len(h))
            mu_ineq = np.zeros(len(g))
        
        return d, lambda_eq, mu_ineq
    
    def _compute_lagrangian_gradient(self, x, grad_f, A_eq, A_ineq, lambda_eq, mu_ineq):
        """ラグランジュ関数の勾配を計算"""
        grad_lag = grad_f(x)
        
        if len(A_eq) > 0 and len(lambda_eq) > 0:
            grad_lag += A_eq.T @ lambda_eq
        
        if len(A_ineq) > 0 and len(mu_ineq) > 0:
            grad_lag += A_ineq.T @ mu_ineq
        
        return grad_lag
    
    def _compute_lagrangian_gradient_difference(self, x1, x2, grad_f, grad_h, grad_g, lambda_eq, mu_ineq):
        """ラグランジュ関数の勾配差を計算（BFGS更新用）"""
        
        # x2でのラグランジュ勾配
        grad_lag_2 = grad_f(x2)
        if grad_h:
            A_eq_2 = np.array([grad_hj(x2) for grad_hj in grad_h])
            if len(lambda_eq) > 0:
                grad_lag_2 += A_eq_2.T @ lambda_eq
        
        if grad_g:
            A_ineq_2 = np.array([grad_gi(x2) for grad_gi in grad_g])
            if len(mu_ineq) > 0:
                grad_lag_2 += A_ineq_2.T @ mu_ineq
        
        # x1でのラグランジュ勾配
        grad_lag_1 = grad_f(x1)
        if grad_h:
            A_eq_1 = np.array([grad_hj(x1) for grad_hj in grad_h])
            if len(lambda_eq) > 0:
                grad_lag_1 += A_eq_1.T @ lambda_eq
        
        if grad_g:
            A_ineq_1 = np.array([grad_gi(x1) for grad_gi in grad_g])
            if len(mu_ineq) > 0:
                grad_lag_1 += A_ineq_1.T @ mu_ineq
        
        return grad_lag_2 - grad_lag_1
    
    def _compute_merit_function(self, x, f, h_vals, g_vals, lambda_eq, mu_ineq, sigma):
        """メリット関数の値を計算"""
        
        f_val = f(x)
        
        if self.merit_function == 'l1':
            penalty = 0
            if len(h_vals) > 0:
                penalty += np.sum(np.abs(h_vals))
            if len(g_vals) > 0:
                penalty += np.sum(np.maximum(0, g_vals))
            return f_val + sigma * penalty
        
        elif self.merit_function == 'l2':
            penalty = 0
            if len(h_vals) > 0:
                penalty += np.sum(h_vals**2)
            if len(g_vals) > 0:
                penalty += np.sum(np.maximum(0, g_vals)**2)
            return f_val + sigma * penalty
        
        elif self.merit_function == 'augmented_lagrangian':
            aug_lag = f_val
            if len(h_vals) > 0 and len(lambda_eq) > 0:
                aug_lag += np.sum(lambda_eq * h_vals)
                aug_lag += (sigma/2) * np.sum(h_vals**2)
            if len(g_vals) > 0 and len(mu_ineq) > 0:
                aug_lag += np.sum(mu_ineq * np.maximum(0, g_vals))
                aug_lag += (sigma/2) * np.sum(np.maximum(0, g_vals)**2)
            return aug_lag
        
        else:
            raise ValueError(f"未知のメリット関数: {self.merit_function}")
    
    def _line_search_merit(self, x, d, f, h, g, lambda_eq, mu_ineq, sigma, 
                          alpha_init=1.0, c1=1e-4, rho=0.5, max_iter=20):
        """メリット関数に基づく直線探索"""
        
        alpha = alpha_init
        
        # 現在のメリット関数値
        h_vals = np.array([hj(x) for hj in h]) if h else np.array([])
        g_vals = np.array([gi(x) for gi in g]) if g else np.array([])
        merit_current = self._compute_merit_function(
            x, f, h_vals, g_vals, lambda_eq, mu_ineq, sigma)
        
        for i in range(max_iter):
            x_trial = x + alpha * d
            
            # 新しい点でのメリット関数値
            h_vals_trial = np.array([hj(x_trial) for hj in h]) if h else np.array([])
            g_vals_trial = np.array([gi(x_trial) for gi in g]) if g else np.array([])
            merit_trial = self._compute_merit_function(
                x_trial, f, h_vals_trial, g_vals_trial, lambda_eq, mu_ineq, sigma)
            
            # 十分減少条件（簡略化版）
            if merit_trial < merit_current:
                return alpha
            
            alpha *= rho
        
        return alpha  # 最小ステップサイズを返す
    
    def _update_hessian_approximation(self, B, s, y):
        """ヘッセ行列近似の更新（BFGS）"""
        
        if self.hessian_update == 'bfgs':
            sy = np.dot(s, y)
            if sy > 1e-12:
                Bs = B @ s
                sBs = np.dot(s, Bs)
                
                B_new = B - np.outer(Bs, Bs) / sBs + np.outer(y, y) / sy
                return B_new
        
        return B

def example_equality_constrained():
    """等式制約付き最適化の例"""
    
    # 問題: min f(x) = (x1-1)² + (x2-2)²
    # subject to: h(x) = x1² + x2² - 1 = 0
    
    def f(x):
        return (x[0] - 1)**2 + (x[1] - 2)**2
    
    def grad_f(x):
        return np.array([2*(x[0] - 1), 2*(x[1] - 2)])
    
    def h1(x):
        return x[0]**2 + x[1]**2 - 1
    
    def grad_h1(x):
        return np.array([2*x[0], 2*x[1]])
    
    h = [h1]
    grad_h = [grad_h1]
    
    # 実行可能な初期点
    x0 = np.array([1.0, 0.0])
    
    # SQP法を実行
    sqp = SQPOptimizer(merit_function='l1')
    x_opt, lambda_opt, mu_opt, history = sqp.optimize(
        f, grad_f, h=h, grad_h=grad_h, x0=x0, verbose=True)
    
    print(f"\n最適解: {x_opt}")
    print(f"最適値: {f(x_opt):.6f}")
    print(f"制約値: {h1(x_opt):.6f}")
    print(f"ラグランジュ乗数: {lambda_opt}")
    
    # 解析解との比較
    # 制約円上での最小点
    direction = np.array([1, 2]) - np.array([0, 0])
    direction_normalized = direction / np.linalg.norm(direction)
    x_analytical = direction_normalized
    
    print(f"解析解: {x_analytical}")
    print(f"誤差: {np.linalg.norm(x_opt - x_analytical):.6f}")
    
    # 結果の可視化
    plt.figure(figsize=(15, 5))
    
    # 収束経路
    plt.subplot(1, 3, 1)
    x_history = np.array(history['x_values'])
    
    # 等高線の描画
    x1_range = np.linspace(-1.5, 1.5, 100)
    x2_range = np.linspace(-1.5, 2.5, 100)
    X1, X2 = np.meshgrid(x1_range, x2_range)
    Z = np.zeros_like(X1)
    
    for i in range(len(x1_range)):
        for j in range(len(x2_range)):
            Z[i,j] = f([X1[i,j], X2[i,j]])
    
    plt.contour(X1, X2, Z, levels=20, alpha=0.6)
    
    # 制約円の描画
    theta = np.linspace(0, 2*np.pi, 100)
    circle_x = np.cos(theta)
    circle_y = np.sin(theta)
    plt.plot(circle_x, circle_y, 'r-', linewidth=2, label='制約: x₁² + x₂² = 1')
    
    # 収束経路
    plt.plot(x_history[:, 0], x_history[:, 1], 'bo-', linewidth=2, markersize=6)
    plt.plot(x_history[0, 0], x_history[0, 1], 'go', markersize=10, label='開始点')
    plt.plot(x_history[-1, 0], x_history[-1, 1], 'ro', markersize=10, label='SQP解')
    plt.plot(x_analytical[0], x_analytical[1], 'k*', markersize=15, label='解析解')
    
    plt.xlabel('x₁')
    plt.ylabel('x₂')
    plt.title('SQP法の収束経路')
    plt.legend()
    plt.grid(True)
    plt.axis('equal')
    
    # 目的関数値の収束
    plt.subplot(1, 3, 2)
    f_values = history['f_values']
    plt.plot(f_values, 'b-o')
    plt.axhline(y=f(x_analytical), color='k', linestyle='--', label='解析解の値')
    plt.xlabel('反復回数')
    plt.ylabel('目的関数値')
    plt.title('目的関数値の収束')
    plt.legend()
    plt.grid(True)
    
    # KKT残差の収束
    plt.subplot(1, 3, 3)
    plt.semilogy(history['kkt_residuals'], 'g-o', label='KKT残差')
    plt.semilogy(history['constraint_violations'], 'r-s', label='制約違反')
    plt.xlabel('反復回数')
    plt.ylabel('残差 (対数スケール)')
    plt.title('収束診断')
    plt.legend()
    plt.grid(True)
    
    plt.tight_layout()
    plt.show()
    
    return x_opt, history

def example_mixed_constraints():
    """等式・不等式混合制約の例"""
    
    # 問題: min f(x) = (x1-2)² + (x2-1)²
    # subject to: h(x) = x1 + x2 - 2 = 0
    #            g(x) = -x1 ≤ 0 (i.e., x1 ≥ 0)
    
    def f(x):
        return (x[0] - 2)**2 + (x[1] - 1)**2
    
    def grad_f(x):
        return np.array([2*(x[0] - 2), 2*(x[1] - 1)])
    
    def h1(x):
        return x[0] + x[1] - 2
    
    def grad_h1(x):
        return np.array([1, 1])
    
    def g1(x):
        return -x[0]  # x1 ≥ 0
    
    def grad_g1(x):
        return np.array([-1, 0])
    
    h = [h1]
    grad_h = [grad_h1]
    g = [g1]
    grad_g = [grad_g1]
    
    # 実行可能な初期点
    x0 = np.array([1.0, 1.0])
    
    sqp = SQPOptimizer(merit_function='l1')
    x_opt, lambda_opt, mu_opt, history = sqp.optimize(
        f, grad_f, h=h, grad_h=grad_h, g=g, grad_g=grad_g, x0=x0, verbose=True)
    
    print(f"\n混合制約問題の最適解: {x_opt}")
    print(f"最適値: {f(x_opt):.6f}")
    print(f"等式制約値: {h1(x_opt):.6f}")
    print(f"不等式制約値: {g1(x_opt):.6f}")
    print(f"等式ラグランジュ乗数: {lambda_opt}")
    print(f"不等式ラグランジュ乗数: {mu_opt}")
    
    return x_opt, history

# 実行例
if __name__ == "__main__":
    print("=== 等式制約付き最適化の例 ===")
    x_opt1, history1 = example_equality_constrained()
    
    print("\n=== 混合制約問題の例 ===")
    x_opt2, history2 = example_mixed_constraints()

MATLAB実装

function [x, lambda, mu, history] = sqp_method(f, grad_f, h, grad_h, g, grad_g, x0, options)
    % SQP法の実装
    
    if nargin < 8, options = struct(); end
    
    % デフォルトパラメータ
    max_iter = getfield_default(options, 'max_iter', 100);
    tol = getfield_default(options, 'tol', 1e-6);
    sigma = getfield_default(options, 'sigma', 1.0);
    
    if isempty(h), h = {}; grad_h = {}; end
    if isempty(g), g = {}; grad_g = {}; end
    
    n = length(x0);
    p = length(h);
    m = length(g);
    
    % 初期化
    x = x0;
    lambda = zeros(p, 1);
    mu = zeros(m, 1);
    B = eye(n);  % 初期ヘッセ行列近似
    
    history.x_values = x0;
    history.f_values = f(x0);
    history.kkt_residuals = [];
    history.constraint_violations = [];
    
    fprintf('SQP法開始: f(x0) = %.6f\n', f(x));
    
    for k = 1:max_iter
        % 制約値の計算
        h_vals = zeros(p, 1);
        for i = 1:p
            h_vals(i) = h{i}(x);
        end
        
        g_vals = zeros(m, 1);
        for i = 1:m
            g_vals(i) = g{i}(x);
        end
        
        % 制約のヤコビアン
        A_eq = zeros(p, n);
        for i = 1:p
            A_eq(i, :) = grad_h{i}(x)';
        end
        
        A_ineq = zeros(m, n);
        for i = 1:m
            A_ineq(i, :) = grad_g{i}(x)';
        end
        
        % 制約違反の計算
        eq_violation = norm(h_vals);
        ineq_violation = sum(max(0, g_vals));
        total_violation = eq_violation + ineq_violation;
        
        % KKT残差の計算
        grad_lag = grad_f(x);
        if p > 0
            grad_lag = grad_lag + A_eq' * lambda;
        end
        if m > 0
            grad_lag = grad_lag + A_ineq' * mu;
        end
        kkt_residual = norm(grad_lag);
        
        % 履歴の更新
        history.kkt_residuals = [history.kkt_residuals, kkt_residual];
        history.constraint_violations = [history.constraint_violations, total_violation];
        
        fprintf('反復 %d: f = %.6f, 制約違反 = %.6f, KKT残差 = %.6f\n', ...
                k, f(x), total_violation, kkt_residual);
        
        % 収束判定
        if kkt_residual < tol && total_violation < tol
            fprintf('\n収束しました (反復回数: %d)\n', k);
            break;
        end
        
        % SQP部分問題を解く
        [d, lambda_new, mu_new] = solve_qp_subproblem(x, grad_f(x), B, ...
                                                     h, grad_h, g, grad_g);
        
        % メリット関数に基づく直線探索
        alpha = line_search_merit(x, d, f, h, g, sigma);
        
        % 更新
        x_new = x + alpha * d;
        
        % ラグランジュ乗数の更新
        if p > 0
            lambda = lambda_new;
        end
        if m > 0
            mu = max(0, mu_new);  % 非負性の確保
        end
        
        % ヘッセ行列近似の更新（BFGS）
        if alpha > 1e-8
            s = x_new - x;
            y = compute_lagrangian_gradient_diff(x, x_new, grad_f, grad_h, grad_g, lambda, mu);
            
            if s' * y > 1e-12
                B = bfgs_update(B, s, y);
            end
        end
        
        x = x_new;
        
        % 履歴の更新
        history.x_values = [history.x_values, x];
        history.f_values = [history.f_values, f(x)];
    end
end

function [d, lambda_new, mu_new] = solve_qp_subproblem(x, grad_fx, B, h, grad_h, g, grad_g)
    % SQP部分問題を解く（簡略化版）
    
    n = length(x);
    p = length(h);
    m = length(g);
    
    % 制約値
    h_vals = zeros(p, 1);
    for i = 1:p
        h_vals(i) = h{i}(x);
    end
    
    g_vals = zeros(m, 1);
    for i = 1:m
        g_vals(i) = g{i}(x);
    end
    
    % 制約のヤコビアン
    A_eq = zeros(p, n);
    for i = 1:p
        A_eq(i, :) = grad_h{i}(x)';
    end
    
    A_ineq = zeros(m, n);
    for i = 1:m
        A_ineq(i, :) = grad_g{i}(x)';
    end
    
    % QP問題を解く（簡単な近似）
    % この部分は実際の実装では専用のQPソルバーを使用
    
    % KKT方程式を解く（等式制約のみの場合）
    if p > 0 && m == 0
        % [B A_eq'; A_eq 0] [d; lambda] = [-grad_fx; -h_vals]
        KKT_matrix = [B, A_eq'; A_eq, zeros(p, p)];
        rhs = [-grad_fx; -h_vals];
        
        if rank(KKT_matrix) == size(KKT_matrix, 1)
            solution = KKT_matrix \ rhs;
            d = solution(1:n);
            lambda_new = solution(n+1:end);
        else
            d = -grad_fx;  % フォールバック
            lambda_new = zeros(p, 1);
        end
        
        mu_new = zeros(m, 1);
    else
        % 一般的なケースは簡略化
        d = -grad_fx;  % 最急降下方向
        lambda_new = zeros(p, 1);
        mu_new = zeros(m, 1);
    end
end

function alpha = line_search_merit(x, d, f, h, g, sigma)
    % L1メリット関数による直線探索
    
    alpha = 1.0;
    rho = 0.5;
    
    % 現在のメリット関数値
    merit_current = compute_l1_merit(x, f, h, g, sigma);
    
    for i = 1:20
        x_trial = x + alpha * d;
        merit_trial = compute_l1_merit(x_trial, f, h, g, sigma);
        
        if merit_trial < merit_current
            return;
        end
        
        alpha = rho * alpha;
    end
end

function merit = compute_l1_merit(x, f, h, g, sigma)
    % L1メリット関数の計算
    
    merit = f(x);
    
    % 等式制約違反
    for i = 1:length(h)
        merit = merit + sigma * abs(h{i}(x));
    end
    
    % 不等式制約違反
    for i = 1:length(g)
        merit = merit + sigma * max(0, g{i}(x));
    end
end

function y = compute_lagrangian_gradient_diff(x1, x2, grad_f, grad_h, grad_g, lambda, mu)
    % ラグランジュ関数の勾配差を計算
    
    % x2でのラグランジュ勾配
    grad_lag_2 = grad_f(x2);
    for i = 1:length(grad_h)
        grad_lag_2 = grad_lag_2 + lambda(i) * grad_h{i}(x2);
    end
    for i = 1:length(grad_g)
        grad_lag_2 = grad_lag_2 + mu(i) * grad_g{i}(x2);
    end
    
    % x1でのラグランジュ勾配
    grad_lag_1 = grad_f(x1);
    for i = 1:length(grad_h)
        grad_lag_1 = grad_lag_1 + lambda(i) * grad_h{i}(x1);
    end
    for i = 1:length(grad_g)
        grad_lag_1 = grad_lag_1 + mu(i) * grad_g{i}(x1);
    end
    
    y = grad_lag_2 - grad_lag_1;
end

function B_new = bfgs_update(B, s, y)
    % BFGS更新
    
    sy = s' * y;
    Bs = B * s;
    sBs = s' * Bs;
    
    B_new = B - (Bs * Bs') / sBs + (y * y') / sy;
end

function val = getfield_default(s, field, default)
    if isfield(s, field)
        val = s.(field);
    else
        val = default;
    end
end

高度なSQP手法

信頼領域SQP

各反復で信頼領域制約を追加：

$$\begin{align} \min_{d \in \mathbb{R}^n} &\quad \frac{1}{2} d^T B_k d + \nabla f(x_k)^T d \\ \text{subject to} &\quad \nabla h_j(x_k)^T d + h_j(x_k) = 0 \\ &\quad \nabla g_i(x_k)^T d + g_i(x_k) \leq 0 \\ &\quad \|d\| \leq \Delta_k \end{align}$$

フィルター法

メリット関数の代わりに、目的関数値と制約違反の組み合わせでフィルターを構成：

• 支配関係: $(f_1, v_1)$ が $(f_2, v_2)$ を支配 ⟺ $f_1 \leq f_2$ かつ $v_1 \leq v_2$
• フィルター: 支配されない点の集合

IP-SQP（内点SQP）

不等式制約をバリア項で処理：

$$\min_{d} \frac{1}{2} d^T B_k d + \nabla f(x_k)^T d - \mu \sum_i \log(-g_i(x_k) - \nabla g_i(x_k)^T d)$$

実用的な考慮事項

QP部分問題の解法

1. アクティブセット法: 小中規模問題に適用
2. 内点法: 大規模問題に効果的
3. 反復法: 非常に大規模な問題

ヘッセ行列近似

1. BFGS更新: 最も一般的
2. 正確なヘッセ行列: 二次収束を保証
3. 有限差分近似: 勾配情報のみ利用可能な場合

大域化戦略

1. メリット関数: L1, L2, 拡張ラグランジュ
2. フィルター法: パレート最適性の概念
3. 信頼領域: ステップサイズの自動調整

まとめ

SQP法の主な特徴：

1. 二次収束性: 適切な条件下でニュートン法と同等
2. 制約処理: 等式・不等式制約の統一的取り扱い
3. 理論的基盤: KKT条件に基づく確固とした理論
4. 実用性: 工学最適化問題での高い成功率
5. 発展性: 様々な改良版が開発されている

成功要因

• 各反復での二次計画問題の正確な解
• 適切なメリット関数の選択
• 堅牢な直線探索またはフィルター法
• 効率的なヘッセ行列近似

制限事項

• QP部分問題の計算コスト
• 初期点の実行可能性（必ずしも必要ではない）
• 大規模問題でのメモリ使用量

実装のポイント

SQP法を実装する際は、QP部分問題の効率的な解法、適切なメリット関数の選択、そして数値的安定性を保つためのヘッセ行列近似が重要です。特に、制約が多い問題では、アクティブセット戦略の実装が成功の鍵となります。

現代的発展

SQP法は現在でも活発に研究されており、大規模問題に対する並列アルゴリズム、不正確なQP解法との組み合わせ、機械学習問題への応用などが注目されています。

参考文献

1. Nocedal, J., & Wright, S. J. (2006). Numerical Optimization. Springer.
2. Boggs, P. T., & Tolle, J. W. (1995). Sequential quadratic programming. Acta Numerica.
3. Fletcher, R. (2013). Practical Methods of Optimization. John Wiley & Sons.
4. Gill, P. E., Murray, W., & Wright, M. H. (2019). Practical Optimization. SIAM.