SciPy (diucapkan sai pie) adalah paket matematika berbasis numpy yang juga mencakup perpustakaan C dan Fortran. SciPy mengubah sesi Python interaktif Anda menjadi lingkungan ilmu data yang lengkap seperti MATLAB, IDL, Octave, R, atau SciLab.
Pada artikel ini, kita akan melihat teknik dasar pemrograman matematika - menyelesaikan masalah optimasi bersyarat untuk fungsi skalar beberapa variabel menggunakan paket scipy.optimize. Algoritme pengoptimalan tanpa batasan telah dibahas di . Bantuan yang lebih detail dan terkini tentang fungsi scipy selalu dapat diperoleh dengan menggunakan perintah help(), Shift+Tab, atau di .
pengenalan
Antarmuka umum untuk memecahkan masalah optimasi bersyarat dan tidak dibatasi dalam paket scipy.optimize disediakan oleh fungsi minimize(). Namun, diketahui bahwa tidak ada metode universal untuk menyelesaikan semua masalah, sehingga pilihan metode yang memadai, seperti biasa, berada di pundak peneliti.
Algoritme optimasi yang sesuai ditentukan menggunakan argumen fungsi minimize(..., method="").
Untuk optimasi bersyarat dari fungsi beberapa variabel, tersedia implementasi metode berikut:
trust-constr— mencari nilai minimum lokal di wilayah kepercayaan. , ;SLSQP— Pemrograman kuadrat sekuensial dengan batasan, metode Newton untuk menyelesaikan sistem Lagrange. .TNC- Newton Terpotong Terkendala, jumlah iterasi terbatas, baik untuk fungsi nonlinier dengan jumlah variabel independen yang banyak. .L-BFGS-B— metode dari tim Broyden–Fletcher–Goldfarb–Shanno, diimplementasikan dengan pengurangan konsumsi memori karena pemuatan sebagian vektor dari matriks Hessian. , .COBYLA— MARE Constrained Optimization By Linear Approximation, optimasi terbatas dengan pendekatan linier (tanpa perhitungan gradien). .
Bergantung pada metode yang dipilih, kondisi dan batasan untuk menyelesaikan masalah ditetapkan secara berbeda:
- objek kelas
Boundsuntuk metode L-BFGS-B, TNC, SLSQP, trust-constr; - Daftar
(min, max)untuk metode yang sama L-BFGS-B, TNC, SLSQP, trust-constr; - suatu objek atau daftar objek
LinearConstraint,NonlinearConstraintuntuk metode COBYLA, SLSQP, trust-constr; - kamus atau daftar kamus
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}untuk COBYLA, metode SLSQP.
Garis besar artikel:
1) Pertimbangkan penggunaan algoritma optimasi bersyarat di wilayah kepercayaan (method=”trust-constr”) dengan batasan yang ditentukan sebagai objek Bounds, LinearConstraint, NonlinearConstraint ;
2) Pertimbangkan pemrograman sekuensial menggunakan metode kuadrat terkecil (method = "SLSQP") dengan batasan yang ditentukan dalam bentuk kamus {'type', 'fun', 'jac', 'args'};
3) Menganalisis contoh optimasi produk manufaktur menggunakan contoh web studio.
Metode pengoptimalan bersyarat="trust-constr"
Implementasi metode trust-constr berdasarkan untuk masalah dengan kendala berupa persamaan dan seterusnya untuk permasalahan yang mempunyai kendala berupa kesenjangan. Kedua metode ini diimplementasikan oleh algoritma untuk menemukan nilai minimum lokal di wilayah kepercayaan dan sangat cocok untuk masalah skala besar.
Rumusan matematis masalah mencari minimum dalam bentuk umum:
Untuk batasan kesetaraan yang ketat, batas bawah ditetapkan sama dengan batas atas .
Untuk batasan satu arah, batas atas atau batas bawah ditetapkan np.inf dengan tanda yang sesuai.
Misalkan kita perlu mencari fungsi minimum Rosenbrock yang diketahui dari dua variabel:
Dalam hal ini, batasan berikut ditetapkan pada domain definisinya:
Dalam kasus kami, ada solusi unik pada saat itu , yang hanya berlaku pada batasan pertama dan keempat.
Mari kita lihat batasannya dari bawah ke atas dan lihat bagaimana kita bisa menulisnya dalam scipy.
Pembatasan и mari kita definisikan menggunakan objek Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Pembatasan и Mari kita tuliskan dalam bentuk linier:
Mari kita definisikan batasan 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 terakhir batasan nonlinier dalam bentuk matriks:
Kami mendefinisikan matriks Jacobian untuk batasan ini dan kombinasi linier dari matriks Hessian dengan vektor arbitrer :
Sekarang kita dapat mendefinisikan batasan nonlinier sebagai sebuah 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 ukurannya besar, matriks juga dapat dispesifikasikan dalam bentuk renggang:
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)Saat menghitung matriks Hessian membutuhkan banyak usaha, Anda dapat menggunakan kelas . Strategi berikut tersedia: BFGS и SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())Hessian juga dapat dihitung menggunakan perbedaan hingga:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')Matriks Jacobian untuk batasan juga dapat dihitung menggunakan perbedaan hingga. Namun dalam hal ini matriks Hessian tidak dapat dihitung dengan menggunakan beda hingga. Hessian harus didefinisikan sebagai fungsi atau menggunakan kelas HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())Solusi dari masalah optimasi terlihat 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 menghitung Hessian dapat didefinisikan 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 produk dari Hessian dan vektor sembarang 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)Alternatifnya, turunan pertama dan kedua dari fungsi yang sedang dioptimalkan dapat diperkirakan. Misalnya, Hessian dapat didekati menggunakan fungsi tersebut SR1 (perkiraan kuasi-Newtonian). Gradien dapat diperkirakan dengan perbedaan hingga.
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)Metode optimasi bersyarat="SLSQP"
Metode SLSQP dirancang untuk menyelesaikan masalah minimalisasi suatu fungsi berupa:
Где и — kumpulan indeks ekspresi yang menggambarkan pembatasan dalam bentuk persamaan atau ketidaksetaraan. — himpunan batas bawah dan atas untuk domain definisi fungsi.
Batasan linier dan nonlinier dijelaskan 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 dilakukan sebagai 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 Optimasi
Sehubungan dengan transisi ke struktur teknologi kelima, mari kita lihat optimasi produksi menggunakan contoh studio web, yang memberi kita pendapatan kecil namun stabil. Bayangkan diri kita sebagai direktur sebuah dapur yang menghasilkan tiga jenis produk:
- x0 - menjual halaman arahan, mulai 10 tr.
- x1 - situs web perusahaan, mulai 20 tr.
- x2 - toko online, mulai 30 tr.
Tim kerja kami yang ramah terdiri dari empat junior, dua menengah, dan satu senior. Dana waktu kerja bulanan mereka:
- Juni:
4 * 150 = 600 чел * час, - tengah:
2 * 150 = 300 чел * час, - senor:
150 чел * час.
Biarkan junior pertama yang tersedia menghabiskan (0, 1, 2) jam untuk pengembangan dan penerapan satu situs tipe (x10, x20, x30), menengah - (7, 15, 20), senior - (5, 10, 15 ) jam waktu terbaik dalam hidup Anda.
Seperti direktur pada umumnya, kami ingin memaksimalkan keuntungan bulanan. Langkah pertama menuju sukses adalah menuliskan fungsi tujuan value sebagai jumlah pendapatan dari produk yang dihasilkan per bulan:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Ini bukan kesalahan; saat mencari nilai maksimum, fungsi tujuan diminimalkan dengan tanda sebaliknya.
Langkah selanjutnya adalah melarang karyawan kami bekerja berlebihan dan menerapkan pembatasan jam kerja:
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]])
}Batasan formalnya adalah keluaran produk hanya boleh positif:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])Dan terakhir, asumsi yang paling menggembirakan adalah karena harga yang murah dan kualitas yang tinggi, antrean pelanggan yang puas terus mengantre untuk kami. Kami dapat memilih sendiri volume produksi bulanan, berdasarkan penyelesaian masalah pengoptimalan yang dibatasi 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 ke bilangan bulat dan hitung beban bulanan pendayung dengan distribusi produk yang optimal x = (8, 6, 3) :
- Juni:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - tengah:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - senor:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
Kesimpulan: agar direktur mendapatkan hasil maksimal yang layak, optimal untuk membuat 8 halaman arahan, 6 situs berukuran sedang, dan 3 toko per bulan. Dalam hal ini, yang senior harus membajak tanpa mengangkat muka dari mesin, beban yang di tengah kira-kira 2/3, yang junior kurang dari setengahnya.
Kesimpulan
Artikel ini menguraikan teknik dasar untuk bekerja dengan paket tersebut scipy.optimize, digunakan untuk menyelesaikan masalah minimisasi bersyarat. Secara pribadi saya menggunakan scipy semata-mata untuk tujuan akademis, itulah sebabnya contoh yang diberikan bersifat lucu.
Banyak sekali teori dan contoh virtual yang dapat ditemukan, misalnya dalam buku karya I.L. Akulich “Pemrograman Matematika dalam Contoh dan Masalah”. Aplikasi yang lebih hardcore scipy.optimize untuk membangun struktur 3D dari sekumpulan gambar () dapat dilihat di .
Sumber informasi utama adalah mereka yang ingin berkontribusi pada terjemahan bagian ini dan bagian lainnya scipy Selamat Datang di .
Terima kasih untuk partisipasi dalam persiapan publikasi.
Sumber: www.habr.com