.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/optimization/constrained/plot_03_penalty_method.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_api_gallery_optimization_constrained_plot_03_penalty_method.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_api_gallery_optimization_constrained_plot_03_penalty_method.py:


The quadratic penalty method: a sequence of unconstrained problems
========================================================================

Fiacco and McCormick's sequential unconstrained minimization replaces
:math:`\min x^2+y^2` subject to :math:`x+y=1` with the unconstrained
problems :math:`\min x^2+y^2+\mu(x+y-1)^2` for growing :math:`\mu`. Each
penalized minimizer violates the constraint slightly, and the
violation shrinks like :math:`1/\mu` as the minimizers approach the
constrained optimum :math:`(0.5, 0.5)`.

.. GENERATED FROM PYTHON SOURCE LINES 14-24

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from mathematicskit.optimization import PenaltyMethod

    f = lambda z: z[0] ** 2 + z[1] ** 2
    grad_f = lambda z: np.array([2.0 * z[0], 2.0 * z[1]])
    h = lambda z: np.array([z[0] + z[1] - 1.0])
    grad_h = lambda z: np.array([[1.0, 1.0]])








.. GENERATED FROM PYTHON SOURCE LINES 25-30

Solve the penalized problems for :math:`\mu = 1, 10, 100, \ldots`
------------------------------------------------------------------
Running the method with one more outer iteration each time exposes the
minimizer of every subproblem. The exact penalized minimizer is
:math:`x = y = \mu/(1+2\mu)`.

.. GENERATED FROM PYTHON SOURCE LINES 30-42

.. code-block:: Python


    x0 = np.array([0.0, 0.0])
    mus, minimizers = [], []
    for n_outer in range(1, 7):
        result = PenaltyMethod(mu0=1.0, mu_factor=10.0, n_outer=n_outer).minimize(f, grad_f, x0, h=h, grad_h=grad_h)
        mus.append(result.extra["mu_history"][-1])
        minimizers.append(result.x)
    mus, minimizers = np.array(mus), np.array(minimizers)
    violation = np.abs(minimizers.sum(axis=1) - 1.0)
    for mu, x, v in zip(mus, minimizers, violation):
        print(f"mu = {mu:>8.0f}: x = {x.round(6)}, exact {mu / (1 + 2 * mu):.6f}, |h(x)| = {v:.1e}")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    mu =        1: x = [0.333333 0.333333], exact 0.333333, |h(x)| = 3.3e-01
    mu =       10: x = [0.47619 0.47619], exact 0.476190, |h(x)| = 4.8e-02
    mu =      100: x = [0.497512 0.497512], exact 0.497512, |h(x)| = 5.0e-03
    mu =     1000: x = [0.49975 0.49975], exact 0.499750, |h(x)| = 5.0e-04
    mu =    10000: x = [0.499975 0.499975], exact 0.499975, |h(x)| = 5.0e-05
    mu =   100000: x = [0.499998 0.499998], exact 0.499998, |h(x)| = 5.0e-06




.. GENERATED FROM PYTHON SOURCE LINES 43-45

Minimizers approach the constraint; violation falls like 1/mu
-------------------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 45-68

.. code-block:: Python


    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4.5))
    xs = np.linspace(-0.2, 1.0, 200)
    X, Y = np.meshgrid(xs, xs)
    for mu, color in zip((1.0, 100.0), ("tab:orange", "tab:green")):
        ax1.contour(X, Y, X**2 + Y**2 + mu * (X + Y - 1.0) ** 2, levels=8, colors=color, linewidths=0.6, alpha=0.7)
    ax1.plot(xs, 1.0 - xs, color="tab:red", lw=2, label=r"constraint $x + y = 1$")
    ax1.plot(*minimizers.T, "o-", color="tab:blue", ms=5, label=r"penalized minimizers $x(\mu)$")
    ax1.plot(0.5, 0.5, "k*", ms=12, label="constrained optimum")
    ax1.set_xlim(-0.2, 1.0)
    ax1.set_ylim(-0.2, 1.0)
    ax1.set_aspect("equal")
    ax1.set_xlabel("x")
    ax1.set_ylabel("y")
    ax1.legend(loc="lower left", fontsize=8)
    ax1.set_title(r"Penalty contours for $\mu = 1$ (orange), $100$ (green)")
    ax2.loglog(mus, violation, "o-", label=r"$|h(x(\mu))|$")
    ax2.loglog(mus, 1.0 / (1.0 + 2.0 * mus), "k--", lw=1, label=r"$1/(1+2\mu)$")
    ax2.set_xlabel(r"penalty weight $\mu$")
    ax2.set_ylabel("constraint violation")
    ax2.legend()
    ax2.set_title("Sequential unconstrained minimization (Fiacco & McCormick)")
    fig.tight_layout()



.. image-sg:: /api/gallery/optimization/constrained/images/sphx_glr_plot_03_penalty_method_001.png
   :alt: Penalty contours for $\mu = 1$ (orange), $100$ (green), Sequential unconstrained minimization (Fiacco & McCormick)
   :srcset: /api/gallery/optimization/constrained/images/sphx_glr_plot_03_penalty_method_001.png
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 0.147 seconds)


.. _sphx_glr_download_api_gallery_optimization_constrained_plot_03_penalty_method.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_03_penalty_method.ipynb <plot_03_penalty_method.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_03_penalty_method.py <plot_03_penalty_method.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_03_penalty_method.zip <plot_03_penalty_method.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_