Metode Penalti

Metode Penalti

--Catatan kuliah MA3071 20210802

Minimumkan $𝑓(𝐱)$ terhadap kendala $𝐱∈𝑆$
Dengan $𝑓$ merupakan fungsi yang kontinu pada $ℝ^𝑛$ dan $𝑆$ adalah himpunan kendala $𝑆=\{𝑥:h_i(𝐱)= 0,i=1,2,...,m\}$

  • Ide dari Metode penalti adalah mengubah masalah optimisasi berkendala menjadi masalah optimisasi tanpa kendala,
  • Kendala tersebut dikalikan dengan suatu konstanta positif (yang membesar) sehingga berfungsi menjadi suatu penalti pada masalah meminimumkan.

Masalah dimodifikasi menjadi :

Minimumkan $𝑓(𝐱)+𝑐𝑃(𝐱)$
Dengan $𝑐$ suatu konstanta positif dan $𝑃(𝐱)$ adalah fungsi pada $ℝ^𝑛$ yang memenuhi

  • $𝑃(𝐱)$ fungsi kontinu
  • $𝑃(𝐱)≥0$ untuk semua $𝐱∈ℝ^𝑛$
    • Biasanya dipilih fungsi penalti $\displaystyle 𝑃(𝐱)=\frac{1}{2} \sum_{i=1}^{m}[ℎ_i(\mathbf{𝑥})]^2$
  • $𝑃(𝐱)=0$ jika dan hanya jika $𝐱∈𝑆$

Contoh

Minimumkan $$f(x_1,x_2)=x_1^2+x_2^2$$ Terhadap kendala $x_1+x_2-1=0$

=======================================================

Dengan metode penalti, masalah menjadi dimodifikasi menjadi,
Minimumkan : $$\theta(\mathbf{x})=x_1^2+x_2^2+\frac{c}{2}\left(x_1+x_2-1\right)^2$$ Dengan $c$ merupakan konstanta positif yang besar $$\nabla \theta(\mathbf{x})=\begin{bmatrix}2x_1+c\left(x_1+x_2-1\right)\\2x_2+c\left(x_1+x_2-1\right)\end{bmatrix}$$ $$\nabla^2 \theta(\mathbf{x})=\begin{bmatrix}2+c&c\\c&2+c\end{bmatrix}$$ Selesaikan masalah optimisasi nonlinear tanpa kendala ini secara numerik

readme first :
-- Operasi aljabar pada python : https://www.ridwanreza.com/2021/07/aljabar-python.html
-- Metode Newton dan Conjugate Gradien : https://www.ridwanreza.com/2021/07/ma3071-algoritma-newton-conjugate.html

In [1]:
import numpy as np

Selesaikan dengan Metode Newton

In [2]:
#Metode Newton

x=np.array([0,0]) #nilai/tebakan awal
print("x0= ",x)
c=1000000 #coba untuk c=1,c=2, c=10, c=100, c=1000000
nablaf=np.array([2*x[0]+c*(x[0]+x[1]-1)  ,  2*x[1]+c*(x[0]+x[1]-1)]) #menghitung nabla f(x) dengan x awal
nabla2f=np.array([[2+c,c],[c,2+c]]) #kebetulan konstan, tidak bergantung pada x
print("nabla theta(x0)= ",nablaf)
print("norm(nabla theta(x0))= ",np.linalg.norm(nablaf))
print("theta(x0)= ",x[0]**2+x[1]**2+c/2*((x[0]+x[1]-1)**2))
error=0.00001

iter=1
maxiter=100

while  np.linalg.norm(nablaf)>error and iter<maxiter:
  x=x-np.linalg.inv(nabla2f)@nablaf #memperbaharui x
  print("\n","x",iter,"= ",x)
  nablaf=np.array([2*x[0]+c*(x[0]+x[1]-1)  ,  2*x[1]+c*(x[0]+x[1]-1)]) #menghitung nabla f(x) dengan x awal
  nabla2f=np.array([[2+c,c],[c,2+c]]) #kebetulan konstan, tidak bergantung pada x
  print("nabla theta(x",iter,")= ",nablaf)
  print("norm(nabla theta(x",iter,")= ",np.linalg.norm(nablaf))
  print("theta(x",iter,")=",x[0]**2+x[1]**2+c/2*((x[0]+x[1]-1)**2))
  iter=iter+1
x0=  [0 0]
nabla theta(x0)=  [-1000000 -1000000]
norm(nabla theta(x0))=  1414213.562373095
theta(x0)=  500000.0

 x 1 =  [0.4999995 0.4999995]
nabla theta(x 1 )=  [-3.62280384e-05 -3.62279529e-05]
norm(nabla theta(x 1 )=  5.12341228091845e-05
theta(x 1 )= 0.49999950000050064

 x 2 =  [0.4999995 0.4999995]
nabla theta(x 2 )=  [-4.98774355e-11 -4.98774355e-11]
norm(nabla theta(x 2 )=  7.053734576184381e-11
theta(x 2 )= 0.4999995000005

Untuk $c=1$ diperoleh nilai $\mathbf{x}=\begin{bmatrix}\frac{1}{4}\\\frac{1}{4}\end{bmatrix}$, titik ini belum memenuhi kendala

Untuk $c=1000000$ diperoleh nilai $\mathbf{x}=\begin{bmatrix}\frac{1}{2}\\\frac{1}{2}\end{bmatrix}$, dengan $f\left(\begin{bmatrix}\frac{1}{2}\\\frac{1}{2}\end{bmatrix}\right)=0.5$

Selesaikan dengan Metode Conjugate Gradien

In [3]:
x=np.array([0,0]) #nilai/tebakan awal
print("x0= ",x)
c=1000000 #coba untuk c=1,c=2, c=10, c=100, c=1000000
nablaf=np.array([2*x[0]+c*(x[0]+x[1]-1)  ,  2*x[1]+c*(x[0]+x[1]-1)]) #menghitung nabla f(x) dengan x awal
print("norm(nabla f(x0))= ",np.linalg.norm(nablaf))
nabla2f=np.array([[2+c,c],[c,2+c]])
d=-nablaf #menghitung d

iter=1
maxiter=100
error=0.000001

while  np.linalg.norm(nablaf)>error and iter<maxiter:
  alpha=-d@nablaf/(d@nabla2f@d)
  x=x+alpha*d
  print("x",iter,"= ", x)
  nablaf=np.array([2*x[0]+c*(x[0]+x[1]-1)  ,  2*x[1]+c*(x[0]+x[1]-1)])
  nabla2f=np.array([[2+c,c],[c,2+c]])
  print("norm(nabla f(x",iter,"))= ",np.linalg.norm(nablaf)) 
  d=-nablaf+( nabla2f@nablaf@d/(d@nabla2f@d)  )*d
  iter=iter+1
x0=  [0 0]
norm(nabla f(x0))=  1414213.562373095
x 1 =  [0.4999995 0.4999995]
norm(nabla f(x 1 ))=  8.647205711577957e-11

Selesaikan dengan Algoritma Umum

Perhatikan apabila masalah tersebut diselesaikan secara numerik degan algoritma optimisasi umum

In [4]:
#Algoritma umum

x=np.array([0,0]) #nilai/tebakan awal
alpha=0.1  #coba kombinasi (0,0,1) , (0,0,0.1), (0,0,0.5), (1,1,1) , (1,1,0.1), (1,1,0.5)
c=1 #coba untuk c=1,c=2, c=10, c=100, c=1000000
print("x0= ",x)
nablaf=np.array([2*x[0]+c*(x[0]+x[1]-1)  ,  2*x[1]+c*(x[0]+x[1]-1)]) #menghitung nabla f(x) dengan x awal

print("nabla f(x0)= ",nablaf)
print("norm(nabla f(x0))= ",np.linalg.norm(nablaf))
error=0.00001

iter=1
maxiter=300

while  np.linalg.norm(nablaf)>error and iter<maxiter:
  x=x-alpha*nablaf #memperbaharui x
  nablaf=np.array([2*x[0]+c*(x[0]+x[1]-1)  ,  2*x[1]+c*(x[0]+x[1]-1)]) #menghitung nabla f(x) yang baru dengan x baru
  if(iter%10==0):
    print("\n","x",iter,"= ",x)
    print("nabla f(x",iter,")= ",nablaf)
    print("norm(nabla f(xx",iter,")= ",np.linalg.norm(nablaf))
  iter=iter+1

print("\n","x",iter-1,"= ",x)
print("nabla f(x",iter-1,")= ",nablaf)
print("norm(nabla f(xx",iter-1,")= ",np.linalg.norm(nablaf))
x0=  [0 0]
nabla f(x0)=  [-1 -1]
norm(nabla f(x0))=  1.4142135623730951

 x 10 =  [0.24848835 0.24848835]
nabla f(x 10 )=  [-0.00604662 -0.00604662]
norm(nabla f(xx 10 )=  0.00855120861640389

 x 20 =  [0.24999086 0.24999086]
nabla f(x 20 )=  [-3.65615844e-05 -3.65615844e-05]
norm(nabla f(xx 20 )=  5.170588852129734e-05

 x 24 =  [0.24999882 0.24999882]
nabla f(x 24 )=  [-4.73838134e-06 -4.73838134e-06]
norm(nabla f(xx 24 )=  6.701083152474188e-06

Apa makna yang dapat diambil dari percobaan di atas?

Apabila ada masukkan silahkan sampaikan

Terima kasih

Comments

Popular posts from this blog

Fungsi Rekursif

Aljabar Python

Pengolahan Matriks (manual)

Contoh Fungsi : Standar Deviasi