SciPy (disebut sai pie) ialah pakej matematik berasaskan numpy yang turut merangkumi perpustakaan C dan Fortran. SciPy menukar sesi Python interaktif anda menjadi persekitaran sains data yang lengkap seperti MATLAB, IDL, Octave, R atau SciLab.
Dalam artikel ini, kita akan melihat teknik asas pengaturcaraan matematik - menyelesaikan masalah pengoptimuman bersyarat untuk fungsi skalar beberapa pembolehubah menggunakan pakej scipy.optimize. Algoritma pengoptimuman tanpa kekangan telah pun dibincangkan dalam . Bantuan yang lebih terperinci dan terkini mengenai fungsi scipy sentiasa boleh diperoleh menggunakan arahan help(), Shift+Tab atau dalam .
Pengenalan
Antara muka biasa untuk menyelesaikan masalah pengoptimuman bersyarat dan tidak terhad dalam pakej scipy.optimize disediakan oleh fungsi minimize(). Walau bagaimanapun, diketahui bahawa tidak ada kaedah universal untuk menyelesaikan semua masalah, jadi pilihan kaedah yang mencukupi, seperti biasa, terletak di bahu penyelidik.
Algoritma pengoptimuman yang sesuai ditentukan menggunakan hujah fungsi minimize(..., method="").
Untuk pengoptimuman bersyarat bagi fungsi beberapa pembolehubah, pelaksanaan kaedah berikut tersedia:
trust-constr— cari minimum tempatan di wilayah keyakinan. , ;SLSQP— pengaturcaraan kuadratik berjujukan dengan kekangan, kaedah Newtonian untuk menyelesaikan sistem Lagrange. .TNC- Newton Terpenggal Terkekang, bilangan lelaran terhad, bagus untuk fungsi tak linear dengan bilangan pembolehubah bebas yang banyak. .L-BFGS-B— kaedah daripada pasukan Broyden–Fletcher–Goldfarb–Shanno, dilaksanakan dengan pengurangan penggunaan ingatan disebabkan pemuatan separa vektor daripada matriks Hessian. , .COBYLA— Pengoptimuman Terkekang MARE Mengikut Penghampiran Linear, pengoptimuman terhad dengan anggaran linear (tanpa pengiraan kecerunan). .
Bergantung pada kaedah yang dipilih, syarat dan sekatan untuk menyelesaikan masalah ditetapkan secara berbeza:
- objek kelas
Boundsuntuk kaedah L-BFGS-B, TNC, SLSQP, trust-constr; - senarai
(min, max)untuk kaedah yang sama L-BFGS-B, TNC, SLSQP, trust-constr; - objek atau senarai objek
LinearConstraint,NonlinearConstraintuntuk kaedah COBYLA, SLSQP, kepercayaan-constr; - kamus atau senarai kamus
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}untuk kaedah COBYLA, SLSQP.
Pelan artikel:
1) Pertimbangkan penggunaan algoritma pengoptimuman bersyarat dalam wilayah amanah (method=”trust-constr”) dengan kekangan yang dinyatakan sebagai objek Bounds, LinearConstraint, NonlinearConstraint ;
2) Pertimbangkan pengaturcaraan berjujukan menggunakan kaedah kuasa dua terkecil (kaedah = "SLSQP") dengan sekatan yang dinyatakan dalam bentuk kamus {'type', 'fun', 'jac', 'args'};
3) Menganalisis contoh pengoptimuman produk perkilangan menggunakan contoh studio web.
Kaedah pengoptimuman bersyarat="trust-constr"
Pelaksanaan kaedah trust-constr berdasarkan untuk masalah dengan kekangan bentuk persamaan dan seterusnya untuk masalah dengan kekangan dalam bentuk ketidaksamaan. Kedua-dua kaedah dilaksanakan oleh algoritma untuk mencari minimum tempatan di rantau keyakinan dan sangat sesuai untuk masalah berskala besar.
Rumusan matematik masalah mencari minimum dalam bentuk umum:
Untuk kekangan kesaksamaan yang ketat, sempadan bawah ditetapkan sama dengan sempadan atas .
Untuk kekangan sehala, had atas atau bawah ditetapkan np.inf dengan tanda yang sepadan.
Biarkan perlu untuk mencari minimum fungsi Rosenbrock yang diketahui bagi dua pembolehubah:
Dalam kes ini, sekatan berikut ditetapkan pada domain definisinya:
Dalam kes kami, terdapat penyelesaian unik pada titik itu , yang mana hanya sekatan pertama dan keempat yang sah.
Mari kita lalui sekatan dari bawah ke atas dan lihat bagaimana kita boleh menulisnya dalam scipy.
Sekatan и mari kita takrifkannya menggunakan objek Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Sekatan и Mari kita tulis dalam bentuk linear:
Mari kita takrifkan kekangan ini sebagai objek LinearConstraint:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])Dan akhirnya kekangan tak linear dalam bentuk matriks:
Kami mentakrifkan matriks Jacobian untuk kekangan ini dan gabungan linear matriks Hessian dengan vektor arbitrari :
Sekarang kita boleh menentukan kekangan tak linear sebagai objek 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)Jika saiznya besar, matriks juga boleh dinyatakan dalam bentuk jarang:
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)atau sebagai objek 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)Apabila mengira matriks Hessian memerlukan banyak usaha, anda boleh menggunakan kelas . Strategi berikut boleh didapati: BFGS и SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())Hessian juga boleh dikira menggunakan perbezaan terhingga:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')Matriks Jacobian untuk kekangan juga boleh dikira menggunakan perbezaan terhingga. Walau bagaimanapun, dalam kes ini matriks Hessian tidak boleh dikira menggunakan perbezaan terhingga. Hessian mesti ditakrifkan sebagai fungsi atau menggunakan kelas HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())Penyelesaian kepada masalah pengoptimuman kelihatan seperti ini:
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]Jika perlu, fungsi untuk mengira Hessian boleh ditakrifkan menggunakan kelas 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)atau hasil darab Hessian dan vektor arbitrari melalui parameter 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)Sebagai alternatif, derivatif pertama dan kedua bagi fungsi yang dioptimumkan boleh dianggarkan. Sebagai contoh, Hessian boleh dianggarkan menggunakan fungsi tersebut SR1 (penghampiran kuasi-Newtonian). Kecerunan boleh dianggarkan dengan perbezaan terhingga.
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)Kaedah pengoptimuman bersyarat = "SLSQP"
Kaedah SLSQP direka untuk menyelesaikan masalah meminimumkan fungsi dalam bentuk:
Где и — set indeks ungkapan yang menerangkan sekatan dalam bentuk kesamaan atau ketaksamaan. — set sempadan bawah dan atas untuk domain takrifan fungsi.
Kekangan linear dan bukan linear diterangkan dalam bentuk kamus dengan kunci 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])
}Pencarian minimum dijalankan seperti berikut:
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 ]Contoh Pengoptimuman
Sehubungan dengan peralihan kepada struktur teknologi kelima, mari kita lihat pengoptimuman pengeluaran menggunakan contoh studio web, yang membawa kita pendapatan yang kecil tetapi stabil. Mari kita bayangkan diri kita sebagai pengarah sebuah kapal dapur yang menghasilkan tiga jenis produk:
- x0 - menjual halaman pendaratan, dari 10 tr.
- x1 - laman web korporat, dari 20 tr.
- x2 - kedai dalam talian, dari 30 tr.
Pasukan kerja mesra kami termasuk empat junior, dua pertengahan dan seorang senior. Dana masa kerja bulanan mereka:
- Jun:
4 * 150 = 600 чел * час, - pertengahan:
2 * 150 = 300 чел * час, - senor:
150 чел * час.
Biarkan junior pertama yang tersedia menghabiskan (0, 1, 2) jam pada pembangunan dan penggunaan satu tapak jenis (x10, x20, x30), tengah - (7, 15, 20), senior - (5, 10, 15 ) jam masa terbaik dalam hidup anda.
Seperti mana-mana pengarah biasa, kami mahu memaksimumkan keuntungan bulanan. Langkah pertama untuk berjaya ialah menulis fungsi objektif value sebagai jumlah pendapatan daripada produk yang dihasilkan sebulan:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Ini bukan ralat; apabila mencari maksimum, fungsi objektif diminimumkan dengan tanda yang bertentangan.
Langkah seterusnya ialah melarang pekerja kami daripada bekerja berlebihan dan memperkenalkan sekatan pada waktu bekerja:
Apa yang setara:
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]])
}Sekatan rasmi ialah keluaran produk mestilah positif sahaja:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])Dan akhirnya, andaian yang paling cerah ialah kerana harga yang rendah dan kualiti yang tinggi, barisan pelanggan yang berpuas hati sentiasa beratur untuk kami. Kami boleh memilih sendiri volum pengeluaran bulanan, berdasarkan penyelesaian masalah pengoptimuman yang terhad 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]Mari kita bulatkan secara longgar kepada nombor bulat dan hitung beban bulanan pendayung dengan pengedaran produk yang optimum x = (8, 6, 3) :
- Jun:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - pertengahan:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - senor:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
Kesimpulan: agar pengarah menerima maksimum yang sepatutnya, adalah optimum untuk membuat 8 halaman pendaratan, 6 tapak bersaiz sederhana dan 3 kedai setiap bulan. Dalam kes ini, senior mesti membajak tanpa melihat ke atas dari mesin, beban pertengahan akan menjadi lebih kurang 2/3, junior kurang daripada separuh.
Kesimpulan
Artikel itu menggariskan teknik asas untuk bekerja dengan pakej scipy.optimize, digunakan untuk menyelesaikan masalah pengecilan bersyarat. Secara peribadi saya gunakan scipy semata-mata untuk tujuan akademik, sebab itu contoh yang diberikan bersifat komik.
Banyak teori dan contoh maya boleh didapati, contohnya, dalam buku oleh I.L. Akulich "Pengaturcaraan matematik dalam contoh dan masalah." Aplikasi yang lebih tegar scipy.optimize untuk membina struktur 3D daripada set imej () boleh dilihat dalam .
Sumber utama maklumat ialah mereka yang ingin menyumbang kepada terjemahan bahagian ini dan bahagian lain scipy Selamat datang ke .
Terima kasih untuk penyertaan dalam penyediaan penerbitan.
Sumber: www.habr.com