Прямо двойственный метод внутренней точки

import numpy as np
 
# Input data
c = np.array([2, 3, 4]) # Objective function coefficients
A = np.array([[1, 1, 0], [2, 0, 1]]) # Constraint matrix
b = np.array([3, 5]) # Constraint vector
x_0 = np.array([0, 0, 0]) # Initial point 
tau = 0.5 # Step size 
epsilon = 10**(-6) # Precision parameter 
 
 
# Auxiliary variables 
n_iterations = 0 # Iteration counter 
x_k = x_0 # Current point 
 
    
# Main loop of the algorithm    
while True:
 
    n_iterations += 1
 
    # Computing the gradient of the objective function at x_k    
    grad_f_xk = c
 
    # Computing the gradient of the dual problem at x_k    
    grad_g_xk = A.T @ (1/b)
 
    # Computing the step direction d_k    
    d_k = -grad_f_xk + tau * grad_g_xk
 
    # Computing the step size alpha    
    alpha = 1
 
    while True:
 
        x_{n+1} = x_{n} + alpha * d_{n}
 
        if c @ x_{n+1} <= c @ x_{n} + epsilon:   break  
 
        alpha /= 2  
 
    if np.linalg.norm(d_{n}) <= epsilon:   break  
 
    x_{n} = x_{n+1}  
 
        
# Output data        
print('Optimal solution is', x_{n+1})        
print('Number of iterations is', n_iterations)

import cvxpy as cp
 
# Define the primal problem variables
x = cp.Variable(n)
 
# Define the objective function
obj = cp.Minimize(cp.sum_squares(A*x - b))
 
# Define the constraints
constraints = [G*x <= h, x >= 0]
 
# Define the primal problem
primal = cp.Problem(obj, constraints)
 
# Define the dual problem variables
y = cp.Variable(m)
 
# Define the objective function
dual_obj = cp.Maximize(b.T*y)
 
# Define the constraints
dual_constraints = [G.T*y >= cp.sum_squares(x), y >= 0]
 
# Define the dual problem
dual = cp.Problem(dual_obj, dual_constraints)
 
# Define the initial feasible point
x_init = cp.Variable(n)
y_init = cp.Variable(m)
init_constraints = [G*x_init <= h, x_init >= 0, G.T*y_init >= cp.sum_squares(x_init), y_init >= 0]
init_problem = cp.Problem(cp.Minimize(0), init_constraints)
init_problem.solve()
 
# Define the parameter for the barrier method
t = cp.Parameter(nonneg=True)
 
# Define the logarithmic barrier function
log_barrier = cp.sum(cp.log(x)) + cp.sum(cp.log(y)) + cp.sum(cp.log(-G*x + h)) + cp.sum(cp.log(-G.T*y + cp.sum_squares(x)))
 
# Define the composite problem
obj_t = obj + t*log_barrier
constraints_t = constraints + [G.T*y == cp.sum_squares(x)]
composite = cp.Problem(obj_t, constraints_t)
 
# Define the parameter for the step size
sigma = cp.Parameter(nonneg=True)
 
# Define the step size and update rule
step = cp.Variable(n+m)
x_step = x + sigma*step[:n]
y_step = y + sigma*step[n:]
 
# Perform the iterations
while t.value > 1e-3:
    # Solve the composite problem
    composite.solve()
    x_star = x.value
    y_star = y.value
 
    # Compute the Newton direction
    grad_f = cp.vstack([2*A.T*(A*x_star - b), -G.T*y_star + cp.sum_squ

import numpy as np
 
def interior_point(c, A, b, x0, s0, alpha=0.1, beta=0.9, tol=1e-6):
    n = x0.shape[0]
    m = b.shape[0]
    x = x0.copy()
    s = s0.copy()
    for k in range(1000):
        # Calculate gradient of the logarithmic barrier function
        grad_L = c - np.dot(A.T, (1 / s))
        # Calculate Jacobian of the constraints
        J = np.dot(A, np.diag(x)) - np.dot(np.diag(s), np.identity(m))
        # Solve the Newton system to get direction
        dx, ds = np.linalg.solve(J, -grad_L)
        # Calculate step size
        t = 1
        while any(x + t * dx <= 0) or any(s + t * ds <= 0):
            t *= beta
        # Update x and s
        x = x + t * dx
        s = s + t * ds
        # Check termination condition
        if np.linalg.norm(grad_L) < tol:
            break
    return x
 
# Example usage
c = np.array([-1, -2])
A = np.array([[1, 1], [2, 1], [1, 2]])
b = np.array([1, 2, 2])
x0 = np.array([1, 1])
s0 = np.array([1, 1, 1])
x = interior_point(c, A, b, x0, s0)
print("Optimal solution:", x)