sesh7.ipynb (4197B)
1 { 2 "cells": [ 3 { 4 "cell_type": "markdown", 5 "id": "f40de83a", 6 "metadata": {}, 7 "source": [ 8 "# Milutin Popovic" 9 ] 10 }, 11 { 12 "cell_type": "code", 13 "execution_count": 1, 14 "id": "5470e45b", 15 "metadata": {}, 16 "outputs": [], 17 "source": [ 18 "import numpy as np\n", 19 "import numdifftools as nd" 20 ] 21 }, 22 { 23 "cell_type": "code", 24 "execution_count": 2, 25 "id": "f0310e79", 26 "metadata": {}, 27 "outputs": [], 28 "source": [ 29 "def lagrange_newton(f, h_s, x_0, mu_0, eps, kmax):\n", 30 " mu_l = [mu_0]; x_l = [x_0]\n", 31 " n = len(x_0); p = len(mu_0)\n", 32 " x_k = x_0; mu_k = mu_0\n", 33 " x_mu_k = np.concatenate([x_k, mu_k])\n", 34 " \n", 35 " f_hess = nd.Hessian(f)\n", 36 " f_grad = nd.Gradient(f)\n", 37 " h_hess_s = [nd.Hessian(h_i) for h_i in h_s]\n", 38 " h_grad_s = [nd.Gradient(h_i) for h_i in h_s]\n", 39 " \n", 40 " for k in range(kmax):\n", 41 " # Construct\n", 42 " d_L_k = f_grad(x_k) + sum(mu_k[i]*h_grad_s[i](x_k) for i in range(p))\n", 43 " dd_L_k = f_hess(x_k) + sum(mu_k[i]*h_hess_s[i](x_k) for i in range(p))\n", 44 " \n", 45 " if p == 1:\n", 46 " h_k = np.array([h_s[0](x_k)])\n", 47 " nabl_h_k = h_grad_s[0](x_k)\n", 48 " part_2 = np.hstack([nabl_h_k, np.zeros(p)])\n", 49 " else:\n", 50 " h_k = np.hstack([h_s[i](x_k) for i in range(p)])\n", 51 " nabl_h_k = np.vstack([h_grad_s[i](x_k) for i in range(p)])\n", 52 " part_2 = np.hstack([nabl_h_k, np.zeros([p, p])])\n", 53 " \n", 54 " part_1 = np.hstack([dd_L_k , nabl_h_k.reshape(n, p)])\n", 55 " nabl_phi = np.vstack([part_1, part_2])\n", 56 " \n", 57 " phi = np.concatenate([d_L_k, h_k])\n", 58 " \n", 59 " if np.linalg.norm(phi) <= eps:\n", 60 " return x_l, mu_l, k\n", 61 " \n", 62 " d_x_mu_k = np.linalg.solve(nabl_phi, -phi)\n", 63 " x_mu_k = x_mu_k + d_x_mu_k\n", 64 " \n", 65 " \n", 66 " x_l.append(x_mu_k[:n])\n", 67 " mu_l.append(x_mu_k[n:])\n", 68 " k += 1\n", 69 " \n", 70 " return x_l, mu_l, kmax" 71 ] 72 }, 73 { 74 "cell_type": "code", 75 "execution_count": 3, 76 "id": "b00a3f83", 77 "metadata": {}, 78 "outputs": [ 79 { 80 "name": "stdout", 81 "output_type": "stream", 82 "text": [ 83 "x_k: [-353.84615385 92.30769231]\n", 84 "μ_k: [1461.53846154]\n", 85 "k : 200\n" 86 ] 87 } 88 ], 89 "source": [ 90 "f_a = lambda x: 2*x[0]**4 + x[1]**4 + 4*x[0]**2 - x[0]*x[1] + 6*x[1]**2\n", 91 "h_a = lambda x: 2*x[0] - x[1] + 4\n", 92 "x_0 = np.array([0, 0]); mu_0 = np.array([0]); kmax = 200; eps = 10e-3\n", 93 "\n", 94 "x_l, mu_l, k = lagrange_newton(f_a, [h_a], x_0, mu_0, eps, kmax)\n", 95 "print('x_k: ', x_l[-1])\n", 96 "print('μ_k: ', mu_l[-1])\n", 97 "print('k :', k)" 98 ] 99 }, 100 { 101 "cell_type": "code", 102 "execution_count": 4, 103 "id": "c6ca49f2", 104 "metadata": {}, 105 "outputs": [ 106 { 107 "name": "stdout", 108 "output_type": "stream", 109 "text": [ 110 "x_k: [-754.14996676 1591.53061069 886.06125938]\n", 111 "μ_k: [435.38303733 469.48428151]\n", 112 "k : 200\n" 113 ] 114 } 115 ], 116 "source": [ 117 "f_b = lambda x: 1000 - x[0]**2 - 2*x[1]**2 - x[2]**2 - x[0]*x[1] - x[0]*x[2]\n", 118 "h_1 = lambda x: x[0]**2 + x[1]**2 + x[2]**2 - 25\n", 119 "h_2 = lambda x: 8*x[0] + 14*x[1] - 7*x[2] - 56\n", 120 "x_0 = np.array([3, 0.2, 3]); mu_0 = np.array([0, 0]); kmax = 200; eps = 10e-5\n", 121 "\n", 122 "x_l, mu_l, k = lagrange_newton(f_b, [h_1, h_2], x_0, mu_0, eps, kmax)\n", 123 "print('x_k: ', x_l[-1])\n", 124 "print('μ_k: ', mu_l[-1])\n", 125 "print('k :', k)" 126 ] 127 } 128 ], 129 "metadata": { 130 "kernelspec": { 131 "display_name": "Python 3 (ipykernel)", 132 "language": "python", 133 "name": "python3" 134 }, 135 "language_info": { 136 "codemirror_mode": { 137 "name": "ipython", 138 "version": 3 139 }, 140 "file_extension": ".py", 141 "mimetype": "text/x-python", 142 "name": "python", 143 "nbconvert_exporter": "python", 144 "pygments_lexer": "ipython3", 145 "version": "3.11.6" 146 } 147 }, 148 "nbformat": 4, 149 "nbformat_minor": 5 150 }