（pythonライブラリ cvxpyを使った) 最適化問題の計算について

Question

前提・実現したいこと

pythonライブラリのcvxpyを使って最適化問題を解こうとしています．
cvxpyを使うのが初めてなので調べながらやっているのですが，
ネット上に乗っていたコードをそのまま実行してもエラーが出るので，
実行できるようにプログラムを変更しています．

まず初めに
プログラム内
x = Variable(n, T + 1)
u = Variable(m, T)
だと以下のエラーが出たので

PS C:\Users\N\Desktop> python test.py
Simulation start
Traceback (most recent call last):
  File "test.py", line 21, in <module>
    x = Variable(n, T + 1)
  File "C:\Users\N\AppData\Local\Programs\Python\Python37\lib\site-packages\cvxpy\expressions\variable.py", line 75, in __init__
    raise TypeError("Variable name %s must be a string." % name)
TypeError: Variable name 51 must be a string.

x = Variable((n, T + 1))
u = Variable((m, T))
に変更しました．

そして実行すると
プログラム上のprob.constraints += [x[:,T] == 0, x[:,0] == x_0]に対してエラーが出てたのですが，他のサイト等を調べても同様に書かれているのでどこが違うのかわかりません．
prob.constraints += [x[:,T] == 0, x[:,0] == x_0]を記述することで制約条件を2つ追加できるそうなのですが．．．．

載せているソースコードは
x = Variable((n, T + 1))
u = Variable((m, T))
に変更したものを載せています．

発生している問題・エラーメッセージ

PS C:\Users\N\Desktop> python test.py
Simulation start
Traceback (most recent call last):
  File "test.py", line 32, in <module>
    prob.constraints += [x[:, T] == 0, x[:, 0] == x_0]
  File "C:\Users\N\AppData\Local\Programs\Python\Python37\lib\site-packages\cvxpy\expressions\expression.py", line 46, in cast_op
    return binary_op(self, other)
  File "C:\Users\N\AppData\Local\Programs\Python\Python37\lib\site-packages\cvxpy\expressions\expression.py", line 575, in __eq__
    return Equality(self, other)
  File "C:\Users\N\AppData\Local\Programs\Python\Python37\lib\site-packages\cvxpy\constraints\zero.py", line 114, in __init__
    self._expr = lhs - rhs
  File "C:\Users\N\AppData\Local\Programs\Python\Python37\lib\site-packages\cvxpy\expressions\expression.py", line 46, in cast_op
    return binary_op(self, other)
  File "C:\Users\N\AppData\Local\Programs\Python\Python37\lib\site-packages\cvxpy\expressions\expression.py", line 464, in __sub__
    return self + -other
  File "C:\Users\N\AppData\Local\Programs\Python\Python37\lib\site-packages\cvxpy\expressions\expression.py", line 46, in cast_op
    return binary_op(self, other)
  File "C:\Users\N\AppData\Local\Programs\Python\Python37\lib\site-packages\cvxpy\expressions\expression.py", line 452, in __add__
    return cvxtypes.add_expr()([self, other])
  File "C:\Users\N\AppData\Local\Programs\Python\Python37\lib\site-packages\cvxpy\atoms\affine\add_expr.py", line 33, in __init__
    super(AddExpression, self).__init__(*arg_groups)
  File "C:\Users\N\AppData\Local\Programs\Python\Python37\lib\site-packages\cvxpy\atoms\atom.py", line 45, in __init__
    self._shape = self.shape_from_args()
  File "C:\Users\N\AppData\Local\Programs\Python\Python37\lib\site-packages\cvxpy\atoms\affine\add_expr.py", line 41, in shape_from_args
    return u.shape.sum_shapes([arg.shape for arg in self.args])
  File "C:\Users\N\AppData\Local\Programs\Python\Python37\lib\site-packages\cvxpy\utilities\shape.py", line 49, in sum_shapes
    len(shapes)*" %s" % tuple(shapes))
ValueError: Cannot broadcast dimensions  (4,) (4, 1)

該当のソースコード

import time
from cvxpy import *
import numpy as np
import matplotlib.pyplot as plt
print("Simulation start")

np.random.seed(1)
n = 4  # state size
m = 2  # input size
T = 50  # number of horizon

# simulation parameter
alpha = 0.2
beta = 5.0

# Model Parameter
A = np.eye(n) + alpha * np.random.randn(n, n)
B = np.random.randn(n, m)
x_0 = beta * np.random.randn(n, 1)

x = Variable((n, T + 1))
u = Variable((m, T))

states = []
for t in range(T):
    cost = sum_squares(x[:, t + 1]) + sum_squares(u[:, t])
    constr = [x[:, t + 1] == A * x[:, t] + B * u[:, t],
              norm(u[:, t], 'inf') <= 1]
    states.append(Problem(Minimize(cost), constr))
# sums problem objectives and concatenates constraints.
prob = sum(states)
prob.constraints += [x[:, T] == 0, x[:, 0] == x_0]

start = time.time()
result = prob.solve(verbose=True)
elapsed_time = time.time() - start
print ("calc time:{0}".format(elapsed_time) + "[sec]")

if result == float("inf"):
    print("Cannot optimize")
    import sys
    sys.exit()
    #  return

#以下グラフの表示
f = plt.figure()
# Plot (u_t)_1.
ax = f.add_subplot(211)
u1 = np.array(u[0, :].value[0, :])[0].tolist()
u2 = np.array(u[1, :].value[0, :])[0].tolist()
plt.plot(u1, '-r', label="u1")
plt.plot(u2, '-b', label="u2")
plt.ylabel(r"$u_t$", fontsize=16)
plt.yticks(np.linspace(-1.0, 1.0, 3))
plt.legend()
plt.grid(True)

# Plot (u_t)_2.
plt.subplot(2, 1, 2)
x1 = np.array(x[0, :].value[0, :])[0].tolist()
x2 = np.array(x[1, :].value[0, :])[0].tolist()
x3 = np.array(x[2, :].value[0, :])[0].tolist()
x4 = np.array(x[3, :].value[0, :])[0].tolist()
plt.plot(range(T + 1), x1, '-r', label="x1")
plt.plot(range(T + 1), x2, '-b', label="x2")
plt.plot(range(T + 1), x3, '-g', label="x3")
plt.plot(range(T + 1), x4, '-k', label="x4")
plt.yticks([-25, 0, 25])
plt.ylim([-25, 25])
plt.ylabel(r"$x_t$", fontsize=16)
plt.xlabel(r"$t$", fontsize=16)
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

試したこと

ここに問題に対して試したことを記載してください。

補足情報（FW/ツールのバージョンなど）

ここにより詳細な情報を記載してください。

Comments

[退会済みユーザー]

x_0のshapeは(4, 1)、x[:, 0]のshapeは(4,)ですので、
prob.constraints += [x[:, T] == 0, x[:, 0] == x_0.flatten()]
とすればshapeの問題は解決できます。
が、これではAttributeError: can't set attributeエラーが出ます。
制約の加え方がわからない、ということでしょうか？

[NR4200]

プログラム内
constr = [x[:, t + 1] == A * x[:, t] + B * u[:, t],
norm(u[:, t], 'inf') <= 1]

で制約条件を入れるのはわかっているのですが，以下のホームページでは
https://myenigma.hatenablog.com/entry/2016/07/25/214014
後から制約条件を加えることができると書かれていたのですが

Answer 1

記事が書かれたのが2016-07-25、当時のcvxpyは0.4.2
現在2020-01-04現在、最新のcvxpyは1.0.25
バージョンがだいぶ変わって文法や用法が変わったものと思われます。

---

私のコメントにありますように、shapeがおかしいことに気づいた人が既に別の問題でこれ関係のバグ報告をしていました。

結果、

• CVXPY v0.4.10では動く
• 新しいバージョンではcreate後の変更は許可されない

だそうで、サンプルの更新が行われたそうです。
一例として別の問題ではありますが動くサンプルが公開されています。

---

ですので、一番早い解決策は古いcvxpyをインストールするです。
古いバージョンにはバグが含まれている可能性がありますので、新しいバージョンのチュートリアルを参考にされた方が良いかと思います。

Comments

[NR4200]

ありがとうございます