Untitled
unknown
plain_text
9 months ago
12 kB
18
Indexable
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
# ---------------------------------------------------------
# PARAMETRAI
# ---------------------------------------------------------
a, b = 7, 0
X0 = np.array([0.0, 0.0])
X1 = np.array([1.0, 1.0])
Xm = np.array([a / 10.0, b / 10.0])
# ---------------------------------------------------------
# TIKSLO FUNKCIJA
# ---------------------------------------------------------
def f(X):
x1, x2 = X
x3 = 1.0 - x1 - x2
if x1 < 0 or x2 < 0:
return 1e5
return -(x1 * x2 * x3) / 8.0
# ---------------------------------------------------------
# GRADIENTAS
# ---------------------------------------------------------
def gradient_f(X):
x1, x2 = X
x3 = 1.0 - x1 - x2
df_dx1 = -(x2 * (x3 - x1)) / 8.0
df_dx2 = -(x1 * (x3 - x2)) / 8.0
return np.array([df_dx1, df_dx2])
# ---------------------------------------------------------
# GRADIENTINIS NUSILEIDIMAS
# ---------------------------------------------------------
def gradient_descent(X0, gamma=3, eps=1e-6, max_iter=2000): # gama=3 optimaliausias is x1, gama=4 optimaliausias is xm
X = X0.copy().astype(float)
hist = [X.copy()]
for i in range(max_iter):
grad = gradient_f(X)
# sustojimo sąlyga pagal gradiento normą
if np.linalg.norm(grad) < eps:
break
X = X - gamma * grad
hist.append(X.copy())
return X, np.array(hist), f(X), i + 1
# ---------------------------------------------------------
# AUKSINIO PJUVIO PAIEŠKA
# ---------------------------------------------------------
def golden_section_search_1d(func, X, d, a0=0, b0=5, tol=1e-6, max_iter=100): #[0,5] optimaliausias is xm [0,4] optimaliausias is x1
tau = (np.sqrt(5) - 1) / 2
l, r = a0, b0
x1 = r - tau * (r - l)
x2 = l + tau * (r - l)
f1 = func(X + x1 * d)
f2 = func(X + x2 * d)
for _ in range(max_iter):
if r - l < tol:
break
if f2 < f1:
l = x1
x1 = x2
f1 = f2
x2 = l + tau * (r - l)
f2 = func(X + x2 * d)
else:
r = x2
x2 = x1
f2 = f1
x1 = r - tau * (r - l)
f1 = func(X + x1 * d)
return (l + r)
# ---------------------------------------------------------
# GREIČIAUSIAS NUSILEIDIMAS
# ---------------------------------------------------------
def steepest_descent(X0, eps=1e-6, max_iter=200):
X = X0.copy().astype(float)
hist = [X.copy()]
for i in range(max_iter):
grad = gradient_f(X)
if np.linalg.norm(grad) < eps:
break
d = -grad
gamma = golden_section_search_1d(f, X, d)
X_new = X + gamma * d
# Ribojimai: x1 ≥ 0, x2 ≥ 0
if (X_new[0] < 0) or (X_new[1] < 0):
gamma = golden_section_search_1d(f, X, d, 0, 1)
X_new = X + gamma * d
if (X_new[0] < 0) or (X_new[1] < 0):
break
X = X_new
hist.append(X.copy())
return X, np.array(hist), f(X), i + 1
# ---------------------------------------------------------
# SIMPLEKSO SUDARYMAS
# ---------------------------------------------------------
def create_regular_simplex(X0, alpha=0.5):
n = len(X0)
delta1 = (np.sqrt(n + 1) + n - 1) / (n * np.sqrt(2)) * alpha
delta2 = (np.sqrt(n + 1) - 1) / (n * np.sqrt(2)) * alpha
simplex = [X0.copy()]
for i in range(n):
Xn = X0.copy()
for j in range(n):
Xn[j] += delta2 if j == i else delta1
simplex.append(Xn)
return np.array(simplex)
# ---------------------------------------------------------
# NELDER–MEAD DEFORMUOJAMO SIMPLEKSO METODAS
# ---------------------------------------------------------
def nelder_mead(X0, alpha=0.5, gamma=2.0, beta=0.5, nu=-0.5, eps=1e-5, max_iter=200):
n = len(X0)
simplex = create_regular_simplex(X0, alpha)
hist = [simplex.copy()]
for it in range(max_iter):
fvals = np.array([f(s) for s in simplex])
idx = np.argsort(fvals)
simplex = simplex[idx]
fvals = fvals[idx]
# --- sustojimo sąlyga ---
f_converged = np.max(np.abs(fvals - fvals[0])) < eps
edge_lengths = [np.linalg.norm(simplex[i] - simplex[j]) for i in range(len(simplex)) for j in range(i + 1, len(simplex))]
shape_converged = np.max(edge_lengths) < eps
if f_converged and shape_converged:
break
Xh = simplex[-1] # blogiausias taškas
Xg = simplex[-2] # antras pagal dydį
Xl = simplex[0] # geriausias taškas
Xc = np.mean(simplex[:-1], axis=0) # vidurkis visų taškų išskyrus blogiausią
# --- 1. Atspindėjimas ---
Xr = Xc + (Xc - Xh)
fr = f(Xr)
if fvals[0] <= fr < fvals[-2]:
# priimame atspindėtą tašką
simplex[-1] = Xr
elif fr < fvals[0]:
# --- 2. Išplėtimas ---
Xe = Xc + gamma * (Xr - Xc)
fe = f(Xe)
if fe < fr:
simplex[-1] = Xe # išplėtimas sėkmingas
else:
simplex[-1] = Xr # tik atspindys
else:
# --- 3. Suspaudimas ---
if fr < fvals[-1]:
# išorinis suspaudimas
Xc_new = Xc + beta * (Xr - Xc)
else:
# vidinis suspaudimas
Xc_new = Xc + nu * (Xh - Xc)
fc_new = f(Xc_new)
if fc_new < fvals[-1]:
simplex[-1] = Xc_new
else:
# --- 4. Simplekso sumažinimas (shrink) ---
for j in range(1, n + 1):
simplex[j] = Xl + 0.5 * (simplex[j] - Xl)
hist.append(simplex.copy())
fvals = np.array([f(s) for s in simplex])
best_idx = np.argmin(fvals)
return simplex[best_idx], f(simplex[best_idx]), it + 1, hist
# ---------------------------------------------------------
# VYKDYMAS IR REZULTATAI
# ---------------------------------------------------------
points = {'X0': X0, 'X1': X1, 'Xm': Xm}
methods = {'Gradientinis': gradient_descent, 'Greičiausias': steepest_descent, 'Nelder-Mead': nelder_mead}
results = {m: {} for m in methods}
for mname, method in methods.items():
for pname, Xs in points.items():
if mname == 'Nelder-Mead':
Xopt, fopt, iters, hist = method(Xs)
results[mname][pname] = {
'X_opt': Xopt,
'f_opt': fopt,
'iters': iters,
'simplex_history': hist,
'final_simplex': hist[-1]
}
else:
Xopt, traj, fopt, iters = method(Xs)
results[mname][pname] = {
'X_opt': Xopt,
'f_opt': fopt,
'iters': iters,
'traj': traj
}
# ---------------------------------------------------------
# SUVESTINĖ (PERKELTA Į VIRŠŲ IR SUSKIRSTYTA PAGAL METODUS)
# ---------------------------------------------------------
print("\n================ GRADIENTINIS NUSILEIDIMAS ================")
print(f"{'Pradinis':<8} {'f(X*)':>12} {'V (apskaičiuota)':>18} {'Iteracijos':>12} {'x1':>12} {'x2':>12}")
for p in results['Gradientinis']:
res = results['Gradientinis'][p]
fopt = res['f_opt']
Xopt = res['X_opt']
iters = res['iters']
V_sq = -fopt if fopt < 1e9 else 0
V = np.sqrt(V_sq) if V_sq > 0 else 0
print(f"{p:<8} {fopt:12.5f} {V:18.5f} {iters:12d} {Xopt[0]:12.5f} {Xopt[1]:12.5f}")
print("\n================ GREIČIAUSIAS NUSILEIDIMAS ================")
print(f"{'Pradinis':<8} {'f(X*)':>12} {'V (apskaičiuota)':>18} {'Iteracijos':>12} {'x1':>12} {'x2':>12}")
for p in results['Greičiausias']:
res = results['Greičiausias'][p]
fopt = res['f_opt']
Xopt = res['X_opt']
iters = res['iters']
V_sq = -fopt if fopt < 1e9 else 0
V = np.sqrt(V_sq) if V_sq > 0 else 0
print(f"{p:<8} {fopt:12.5f} {V:18.5f} {iters:12d} {Xopt[0]:12.5f} {Xopt[1]:12.5f}")
print("\n================ NELDER–MEAD METODAS ================")
print(f"{'Pradinis':<8} {'f(X*)':>12} {'V (apskaičiuota)':>18} {'Iteracijos':>12} {'x1':>12} {'x2':>12}")
for p in results['Nelder-Mead']:
res = results['Nelder-Mead'][p]
fopt = res['f_opt']
Xopt = res['X_opt']
iters = res['iters']
V_sq = -fopt if fopt < 1e9 else 0
V = np.sqrt(V_sq) if V_sq > 0 else 0
print(f"{p:<8} {fopt:12.5f} {V:18.5f} {iters:12d} {Xopt[0]:12.5f} {Xopt[1]:12.5f}")
# ---------------------------------------------------------
# GRAFIKAI
# ---------------------------------------------------------
x1_range = np.linspace(0.01, 0.98, 120)
x2_range = np.linspace(0.01, 0.98, 120)
X1g, X2g = np.meshgrid(x1_range, x2_range)
Fgrid = np.full_like(X1g, np.nan)
for i in range(X1g.shape[0]):
for j in range(X1g.shape[1]):
x1, x2 = X1g[i, j], X2g[i, j]
if x1 + x2 <= 1 and x1 >= 0 and x2 >= 0:
Fgrid[i, j] = f([x1, x2])
fig, axes = plt.subplots(3, 3, figsize=(15, 15))
axes = axes.flatten()
colors = {'X0': 'red', 'X1': 'blue', 'Xm': 'green'}
# Gradientinis nusileidimas (1–3)
pos_grad = {'X0': 0, 'X1': 1, 'Xm': 2}
for pname, pos in pos_grad.items():
ax = axes[pos]
ax.contour(X1g, X2g, Fgrid, levels=25, cmap='viridis')
traj_gr = results['Gradientinis'][pname]['traj']
ax.plot(traj_gr[:, 0], traj_gr[:, 1], 'o-', color=colors[pname], markersize=3)
ax.set_title(f'Gradientinis - {pname}')
ax.set_xlim(0, 1)
ax.set_ylim(0, 1)
ax.set_aspect('equal')
ax.set_xlabel('x1')
ax.set_ylabel('x2')
# Greičiausias nusileidimas (4–6)
pos_steep = {'X0': 3, 'X1': 4, 'Xm': 5}
for pname, pos in pos_steep.items():
ax = axes[pos]
ax.contour(X1g, X2g, Fgrid, levels=25, cmap='viridis')
traj = results['Greičiausias'][pname]['traj']
ax.plot(traj[:, 0], traj[:, 1], 'o-', color=colors[pname], markersize=3)
ax.set_title(f'Greičiausias - {pname}')
ax.set_xlim(0, 1)
ax.set_ylim(0, 1)
ax.set_aspect('equal')
ax.set_xlabel('x1')
ax.set_ylabel('x2')
# Nelder-Mead (7–9)
plot_positions_nm = {'X0':6, 'X1':7, 'Xm':8}
for pname, pos in plot_positions_nm.items():
ax = axes[pos]
ax.contour(X1g, X2g, Fgrid, levels=25, cmap='viridis')
nm_data = results['Nelder-Mead'][pname]
simplex_history = nm_data['simplex_history']
# Braižome visus simpleksus evoliucijos metu
for k, simp in enumerate(simplex_history):
tri = np.vstack((simp, simp[0]))
ax.plot(tri[:,0], tri[:,1], '-',
linewidth=max(0.5, 1.5*(1 - 0.9*k/len(simplex_history))),
alpha=max(0.15, 1 - 0.85*k/len(simplex_history)), color='orange')
# Galutinį simpleksą paryškiname raudonai
final_simplex = nm_data['final_simplex']
ax.fill(final_simplex[:,0], final_simplex[:,1], color='red', alpha=0.3, label='galutinis simpleksas')
final_tri = np.vstack((final_simplex, final_simplex[0]))
ax.plot(final_tri[:,0], final_tri[:,1], 'r-', linewidth=2)
# Pradinis simpleksas
initial_simplex = create_regular_simplex(points[pname], alpha=0.5)
ax.plot(initial_simplex[:,0], initial_simplex[:,1], 'x', color='black', label='initial simplex verts')
# Optimalus taškas
best = nm_data['X_opt']
ax.plot(best[0], best[1], 'X', color='magenta', markersize=10, label='opt')
ax.set_xlim(0, 1)
ax.set_ylim(0, 1)
ax.set_aspect('equal')
ax.set_title(f'Nelder-Mead - {pname}')
ax.set_xlabel('x1')
ax.set_ylabel('x2')
ax.legend(fontsize=8)
plt.tight_layout()
plt.show()
Editor is loading...
Leave a Comment