14

CVXOPTでの学習課題として以下のことを試みています。不等式制約を削除し、等式制約をいくつか追加することで、ここのコード例に小さな変更を加えました。

from cvxopt import solvers, blas, matrix, spmatrix, spdiag, log, div
solvers.options['show_progress'] = False
import numpy as np    
np.random.seed(1)

# minimize     p'*log(p)
# subject to
#              sum(p) = 1
#              sum(p'*a) = target1
#              sum(p'*max(a-K,a^2)) = target2

a = np.random.randint(20, 30, size=500)
target1 = 30
target2 = 0.60
K = 26

A = matrix(np.vstack([np.ones(500), a, np.array([max(x-K,x*x) for x in a])]))
b = matrix([1.0, target1, target2])

n = 500
def F(x=None, z=None):
   if x is None: return 0, matrix(1.0, (n,1))
   if min(x) <= 0: return None
   f = x.T*log(x)
   grad = 1.0 + log(x)
   if z is None: return f, grad.T
   H = spdiag(z[0] * x**-1)
   return f, grad.T, H
sol = solvers.cp(F, A=A, b=b)
p = sol['x']

しかし、次のことを実行すると：

np.sum(p)
243.52686763225338

これは、最適化の最初の制約に違反しています。ここで何が問題なのかわかりません。(乱数を使用して変数を生成しているため、異なる値が生成されることに注意してください。ただし、私aのnp.sum(p)ものと同じ違反を観察する必要があります。

元のリンクの不等式制約を保持し、さらに 2 つの等式制約を追加しても、等式制約に違反します。

確実に使用できる他のパッケージ、つまり維持されているパッケージはありますか?

編集:実行可能な解決策がない場合、実行可能な解決策が見つからないというメッセージが表示されるべきではありませんか?

4

1 に答える 1

5

@tihomがコメントしたように、問題は実行不可能です。これが解決したい問題であると本当に確信していますか? 最初の制約は次のことを意味します。

p1 + p2 + ... + pn = 1
p1*a1 + p2*a2 + ... + an*pn = 30
p1*a1^2 + p2*a2^2 + ... pn*an^2 = 0.6

ai >= 20for all であるため、最後の制約が最初または 2 番目の制約と同時に満たされることはありませんi。つまり、合計p1*a1^2 + p2*a2^2 + ... pn*an^2は常に他の合計よりも大きくなります ( に注意してくださいpi > 0)。

制約を満たす点があれば、最適解target1 = sum(a/500.)を見つけることができます。target2= sum(a*a/500.)

最後の制約では、最大値が に単純化されることに注意してください。max(a - K, a^2) = a^2これは、 に関係なく真ですa。

編集: ソリューション (例: print sol) を調べると、次のような結果が得られます。

{'status': 'unknown', 'zl': <0x1 matrix, tc='d'>, 'dual slack': 1.0000000000000007, 'relative gap': 0.005911420508296136, 'dual objective': -97.17320604198335, 'snl': <0x1 matrix, tc='d'>, 'gap': 0.9737154924375709, 'primal objective': -164.7176835197311, 'primal slack': 0.9737154924375703, 'znl': <0x1 matrix, tc='d'>, 'primal infeasibility': 0.5114570271204905, 'dual infeasibility': 0.5091221046374248, 'sl': <0x1 matrix, tc='d'>, 'y': <3x1 matrix, tc='d'>, 'x': <500x1 matrix, tc='d'>}

がであることに注意してください。statusつまり'unknown'、実行可能な解が見つかりません。これは次のドキュメントにありcvxopt.solvers.cpます: http://cvxopt.org/userguide/solvers.html?highlight=cp#cvxopt.solvers.cp

于 2016-12-31T19:36:01.273 に答える