SciPy, 조건에 따른 최적화

SciPy(sai Pie로 발음)는 C 및 Fortran 라이브러리도 포함하는 numpy 기반 수학 패키지입니다. SciPy는 대화형 Python 세션을 MATLAB, IDL, Octave, R 또는 SciLab과 같은 완전한 데이터 과학 환경으로 전환합니다.

이 기사에서는 scipy.optimize 패키지를 사용하여 여러 변수의 스칼라 함수에 대한 조건부 최적화 문제를 해결하는 수학적 프로그래밍의 기본 기술을 살펴보겠습니다. 제약 없는 최적화 알고리즘은 이미 다음에서 논의되었습니다. 마지막 기사. scipy 기능에 대한 더 자세한 최신 도움말은 help() 명령, Shift+Tab 또는 공식 문서.

소개

scipy.optimize 패키지의 조건부 최적화 문제와 제약 없는 최적화 문제를 모두 해결하기 위한 공통 인터페이스는 다음 함수로 제공됩니다. minimize(). 그러나 모든 문제를 해결하는 보편적인 방법은 없다는 것이 알려져 있으므로, 언제나 그렇듯이 적절한 방법을 선택하는 것은 연구자의 몫입니다.
적절한 최적화 알고리즘은 함수 인수를 사용하여 지정됩니다. minimize(..., method="").
여러 변수의 함수에 대한 조건부 최적화를 위해 다음 방법을 구현할 수 있습니다.

• trust-constr — 신뢰 영역에서 지역 최소값을 검색합니다. 위키 기사, 하브레에 관한 기사;
• SLSQP — 제약 조건이 있는 순차 XNUMX차 계획법, 라그랑주 시스템을 해결하기 위한 뉴턴 방식. 위키 기사.
• TNC - Truncated Newton 제약이 있고 반복 횟수가 제한되어 있으며 독립 변수가 많은 비선형 함수에 적합합니다. 위키 기사.
• L-BFGS-B — Broyden–Fletcher–Goldfarb–Shanno 팀의 방법으로 헤세 행렬에서 벡터를 부분적으로 로드하여 메모리 소비를 줄이도록 구현되었습니다. 위키 기사, 하브레에 관한 기사.
• COBYLA — 선형 근사를 통한 MARE 제약 최적화, 선형 근사를 통한 제약 최적화(그라디언트 계산 없음). 위키 기사.

선택한 방법에 따라 문제 해결을 위한 조건과 제한 사항이 다르게 설정됩니다.

• 클래스 개체 Bounds L-BFGS-B, TNC, SLSQP, trust-constr 방법의 경우;
• 목록 (min, max) 동일한 방법으로 L-BFGS-B, TNC, SLSQP, trust-constr;
• 객체 또는 객체 목록 LinearConstraint, NonlinearConstraint COBYLA, SLSQP, trust-constr 방법의 경우;
• 사전 또는 사전 목록 {'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt} COBYLA, SLSQP 방법의 경우.

기사 개요:
1) 객체로 지정된 제약 조건을 사용하여 신뢰 영역(method="trust-constr")에서 조건부 최적화 알고리즘의 사용을 고려합니다. Bounds, LinearConstraint, NonlinearConstraint ;
2) 사전 형식으로 지정된 제한 사항을 사용하여 최소 제곱법(방법 = "SLSQP")을 사용하는 순차 프로그래밍을 고려합니다. {'type', 'fun', 'jac', 'args'};
3) 웹스튜디오를 사례로 제작된 제품의 최적화 사례를 분석한다.

조건부 최적화 방법="trust-constr"

방법의 구현 trust-constr 기반으로 EQSQP 평등의 형태에 대한 제약과 관련된 문제 TRIP 불평등 형태의 제약 조건이 있는 문제의 경우. 두 방법 모두 신뢰 영역에서 국소 최소값을 찾는 알고리즘으로 구현되며 대규모 문제에 매우 적합합니다.

일반적인 형태로 최소값을 찾는 문제의 수학적 공식화:

엄격한 동등 제약 조건의 경우 하한은 상한과 동일하게 설정됩니다. .
단방향 제약조건의 경우 상한 또는 하한이 설정됩니다. np.inf 해당 기호가 있습니다.
알려진 두 변수의 Rosenbrock 함수의 최소값을 찾는 것이 필요하다고 가정합니다.

이 경우 정의 영역에 다음과 같은 제한 사항이 설정됩니다.

우리의 경우에는 그 시점에 독특한 솔루션이 있습니다. , 첫 번째와 네 번째 제한사항만 유효합니다.
제한 사항을 아래에서 위로 살펴보고 scipy에서 이를 작성하는 방법을 살펴보겠습니다.
제한 и Bounds 객체를 사용하여 정의해 보겠습니다.

from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])

제한 и 선형 형식으로 작성해 보겠습니다.

이러한 제약 조건을 LinearConstraint 개체로 정의해 보겠습니다.

import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])

마지막으로 행렬 형식의 비선형 제약 조건은 다음과 같습니다.

이 제약 조건에 대해 야코비안 행렬을 정의하고 헤세 행렬과 임의 벡터의 선형 결합을 정의합니다. :

이제 비선형 제약 조건을 객체로 정의할 수 있습니다. NonlinearConstraint:

from scipy.optimize import NonlinearConstraint

def cons_f(x):
     return [x[0]**2 + x[1], x[0]**2 - x[1]]

def cons_J(x):
     return [[2*x[0], 1], [2*x[0], -1]]

def cons_H(x, v):
     return v[0]*np.array([[2, 0], [0, 0]]) + v[1]*np.array([[2, 0], [0, 0]])

nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=cons_H)

크기가 큰 경우 행렬을 희소 형식으로 지정할 수도 있습니다.

from scipy.sparse import csc_matrix

def cons_H_sparse(x, v):
     return v[0]*csc_matrix([[2, 0], [0, 0]]) + v[1]*csc_matrix([[2, 0], [0, 0]])

nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1,
                                            jac=cons_J, hess=cons_H_sparse)

아니면 객체로 LinearOperator:

from scipy.sparse.linalg import LinearOperator

def cons_H_linear_operator(x, v):
    def matvec(p):
        return np.array([p[0]*2*(v[0]+v[1]), 0])
    return LinearOperator((2, 2), matvec=matvec)

nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1,
                                jac=cons_J, hess=cons_H_linear_operator)

헤세 행렬을 계산할 때 많은 노력이 필요하므로 수업을 사용할 수 있습니다 HessianUpdateStrategy. 다음 전략을 사용할 수 있습니다. BFGS и SR1.

from scipy.optimize import BFGS

nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())

헤세 행렬은 유한 차분을 사용하여 계산할 수도 있습니다.

nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')

제약 조건에 대한 야코비 행렬은 유한 차분을 사용하여 계산할 수도 있습니다. 그러나 이 경우 헤세 행렬은 유한 차분을 사용하여 계산할 수 없습니다. Hessian은 함수로 정의되거나 HessianUpdateStrategy 클래스를 사용하여 정의되어야 합니다.

nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())

최적화 문제에 대한 해결책은 다음과 같습니다.

from scipy.optimize import minimize
from scipy.optimize import rosen, rosen_der, rosen_hess, rosen_hess_prod

x0 = np.array([0.5, 0])
res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hess=rosen_hess,
                constraints=[linear_constraint, nonlinear_constraint],
                options={'verbose': 1}, bounds=bounds)
print(res.x)

`gtol` termination condition is satisfied.
Number of iterations: 12, function evaluations: 8, CG iterations: 7, optimality: 2.99e-09, constraint violation: 1.11e-16, execution time: 0.033 s.
[0.41494531 0.17010937]

필요한 경우 헤세 행렬을 계산하는 함수는 LinearOperator 클래스를 사용하여 정의할 수 있습니다.

def rosen_hess_linop(x):
    def matvec(p):
        return rosen_hess_prod(x, p)
    return LinearOperator((2, 2), matvec=matvec)

res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hess=rosen_hess_linop,
                 constraints=[linear_constraint, nonlinear_constraint],
                 options={'verbose': 1}, bounds=bounds)

print(res.x)

또는 헤세 행렬과 매개변수를 통한 임의의 벡터의 곱 hessp:

res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hessp=rosen_hess_prod,
                constraints=[linear_constraint, nonlinear_constraint],
                options={'verbose': 1}, bounds=bounds)
print(res.x)

대안적으로, 최적화되는 함수의 XNUMX차 도함수와 XNUMX차 도함수를 근사화할 수 있습니다. 예를 들어, 헤세 행렬은 다음 함수를 사용하여 근사화할 수 있습니다. SR1 (준뉴턴 근사). 그래디언트는 유한 차분으로 근사화될 수 있습니다.

from scipy.optimize import SR1
res = minimize(rosen, x0, method='trust-constr',  jac="2-point", hess=SR1(),
               constraints=[linear_constraint, nonlinear_constraint],
               options={'verbose': 1}, bounds=bounds)
print(res.x)

조건부 최적화 방법="SLSQP"

SLSQP 방법은 다음 형식으로 함수를 최소화하는 문제를 해결하도록 설계되었습니다.

Где и — 동등 또는 부등의 형태로 제한을 설명하는 표현의 인덱스 세트. — 함수 정의 영역의 하한 및 상한 세트.

선형 및 비선형 제약 조건은 키가 있는 사전 형식으로 설명됩니다. type, fun и jac.

ineq_cons = {'type': 'ineq',
             'fun': lambda x: np.array ([1 - x [0] - 2 * x [1],
                                          1 - x [0] ** 2 - x [1],
                                          1 - x [0] ** 2 + x [1]]),
             'jac': lambda x: np.array ([[- 1.0, -2.0],
                                          [-2 * x [0], -1.0],
                                          [-2 * x [0], 1.0]])
            }

eq_cons = {'type': 'eq',
           'fun': lambda x: np.array ([2 * x [0] + x [1] - 1]),
           'jac': lambda x: np.array ([2.0, 1.0])
          }

최소값 검색은 다음과 같이 수행됩니다.

x0 = np.array([0.5, 0])
res = minimize(rosen, x0, method='SLSQP', jac=rosen_der,
               constraints=[eq_cons, ineq_cons], options={'ftol': 1e-9, 'disp': True},
               bounds=bounds)

print(res.x)

Optimization terminated successfully.    (Exit mode 0)
            Current function value: 0.34271757499419825
            Iterations: 4
            Function evaluations: 5
            Gradient evaluations: 4
[0.41494475 0.1701105 ]

최적화 예시

XNUMX차 기술구조로의 전환과 관련하여 작지만 안정적인 수입을 가져다주는 웹스튜디오를 예로 들어 생산 최적화를 살펴보자. 우리 자신이 세 가지 유형의 제품을 생산하는 주방의 책임자라고 상상해 봅시다.

• x0 - 랜딩 페이지 판매, 10 tr부터.
• x1 - 기업 웹사이트, 20 tr.
• x2 - 30 tr부터 온라인 상점.

우리의 친절한 작업 팀에는 후배 XNUMX명, 중간 XNUMX명, 선배 XNUMX명이 포함되어 있습니다. 월간 근무 시간 기금:

• XNUMX월: 4 * 150 = 600 чел * час,
• 중간: 2 * 150 = 300 чел * час,
• 선생님: 150 чел * час.

사용 가능한 첫 번째 주니어가 (x0, x1, x2), 중간 - (10, 20, 30), 시니어 - (7, 15, 20) 유형의 사이트 하나를 개발 및 배포하는 데 (5, 10, 15) 시간을 소비하도록 합니다. ) 인생에서 가장 좋은 시간의 시간.

여느 이사와 마찬가지로 우리도 월 수익을 극대화하고 싶습니다. 성공의 첫 번째 단계는 목적 함수를 작성하는 것입니다. value 한 달에 생산된 제품의 수입 금액:

def value(x):
    return - 10*x[0] - 20*x[1] - 30*x[2]

이는 오류가 아닙니다. 최대값을 검색할 때 목적 함수는 반대 부호로 최소화됩니다.

다음 단계는 직원의 과로를 금지하고 근무 시간에 대한 제한을 도입하는 것입니다.

동등한 것:

ineq_cons = {'type': 'ineq',
             'fun': lambda x: np.array ([600 - 10 * x [0] - 20 * x [1] - 30 * x[2],
                                         300 - 7  * x [0] - 15 * x [1] - 20 * x[2],
                                         150 - 5  * x [0] - 10 * x [1] - 15 * x[2]])
            }

공식적인 제한은 제품 출력이 긍정적이어야 한다는 것입니다.

bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])

그리고 마지막으로, 가장 장밋빛 가정은 저렴한 가격과 높은 품질로 인해 만족한 고객들이 끊임없이 우리를 기다리고 있다는 것입니다. 제한된 최적화 문제를 해결하여 월간 생산량을 직접 선택할 수 있습니다. scipy.optimize:

x0 = np.array([10, 10, 10])
res = minimize(value, x0, method='SLSQP', constraints=ineq_cons, bounds=bnds)
print(res.x)

[7.85714286 5.71428571 3.57142857]

정수로 느슨하게 반올림하고 제품의 최적 분포를 통해 조정자의 월간 부하를 계산해 보겠습니다. x = (8, 6, 3) :

• XNUMX월: 8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час;
• 중간: 8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час;
• 선생님: 8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.

결론: 이사가 마땅한 최대치를 얻으려면 매월 랜딩 페이지 8개, 중간 규모 사이트 6개, 매장 3개를 만드는 것이 최적입니다. 이 경우 선배는 기계에서 눈을 떼지 않고 쟁기질을 해야 하며, 중간의 부하는 약 2/3, 후배는 절반 미만이 됩니다.

결론

이 기사에서는 패키지 작업을 위한 기본 기술을 설명합니다. scipy.optimize, 조건부 최소화 문제를 해결하는 데 사용됩니다. 개인적으로 나는 사용한다 scipy 순전히 학문적 목적으로 제공된 예가 그토록 코믹한 이유입니다.

예를 들어, I.L. Akulich의 책 "예제와 문제의 수학 프로그래밍"에서 많은 이론과 가상 예를 찾을 수 있습니다. 더 많은 하드코어 애플리케이션 scipy.optimize 이미지 세트에서 3D 구조를 구축합니다(하브레에 관한 기사)에서 볼 수 있습니다. scipy 요리책.

주요 정보 출처는 docs.scipy.org이 섹션과 다른 섹션의 번역에 기여하고 싶은 사람들 scipy 에 오신 것을 환영합니다 GitHub의.

감사 메피스토피 출판물 준비에 참여하기 위해.

출처 : habr.com