静的解析のための for ループの制約の変換

ライブ スクリプトを開く

最適化式の静的解析で説明されているように、静的解析は問題ベースの式の for ループを高速化します。静的解析は、比較演算子 <=、==、>= ではなく、fcn2optimexpr によって作成された式に適用されます。そのため、for ループを使用して作成した制約がある場合は、比較の両辺から制約の右辺を引き、比較がゼロ、あるいは定数になるようにします。次に、fcn2optimexpr を適用します。

この例では、i=1,2,...,20 に対して xi≤u+i-1 と xi≥-u-i+1 の制約があります。まず最適化変数を作成します。

N = 20;
x = optimvar("x",N);
u = optimvar("u");

for ループの制約を表すには、制約が 0 と比較されるように適切な値を引きます。

xi-u-i+1≤0xi+u+i-1≥0.

通常、これらの制約は次のコードで表します。

for i = 1:N
    cons1(i) = x(i) - u - i + 1;
    cons2(i) = x(i) + u + i - 1;
end

静的解析を適用するには、この for ループを forloopfcn という名前の別の補助関数に配置します。結果として得られる関数は、この例の終わりに掲載しています。(静的解析のための for ループの作成に示すように、スクリプト内の for ループを別個の関数に容易に変換できます。)

静的解析を活用するには、fcn2optimexpr を使用して関数を最適化式に変換します。

[cons1,cons2] = fcn2optimexpr(@forloopfcn,x,u);

最適化問題を作成し、制約を問題に配置します。

prob = optimproblem;
prob.Constraints.cons1 = cons1 <= 0;
prob.Constraints.cons2 = cons2 >= 0;

変数 u の重み付き和と、ベクトル 1.5*(1:N) と変数 x の差の二乗和である目的関数を作成します。

M = 100; % Weight factor for u
prob.Objective = M*u + sum((x - 1.5*(1:N)').^2);

問題を確認します。

show(prob)
  OptimizationProblem : 

	Solve for:
       u, x

	minimize :
       x(1)^2 + x(2)^2 + x(3)^2 + x(4)^2 + x(5)^2 + x(6)^2 + x(7)^2 + x(8)^2 + x(9)^2 + x(10)^2 + x(11)^2 + x(12)^2 + x(13)^2 + x(14)^2 + x(15)^2 + x(16)^2 + x(17)^2 + x(18)^2 + x(19)^2 + x(20)^2 + 100*u - 3*x(1) - 6*x(2) - 9*x(3) - 12*x(4) - 15*x(5) - 18*x(6) - 21*x(7) - 24*x(8) - 27*x(9) - 30*x(10) - 33*x(11) - 36*x(12) - 39*x(13) - 42*x(14) - 45*x(15) - 48*x(16) - 51*x(17) - 54*x(18) - 57*x(19) - 60*x(20) + 6457.5


	subject to cons1:
       -u + x(1) <= 0
       -u + x(2) <= 1
       -u + x(3) <= 2
       -u + x(4) <= 3
       -u + x(5) <= 4
       -u + x(6) <= 5
       -u + x(7) <= 6
       -u + x(8) <= 7
       -u + x(9) <= 8
       -u + x(10) <= 9
       -u + x(11) <= 10
       -u + x(12) <= 11
       -u + x(13) <= 12
       -u + x(14) <= 13
       -u + x(15) <= 14
       -u + x(16) <= 15
       -u + x(17) <= 16
       -u + x(18) <= 17
       -u + x(19) <= 18
       -u + x(20) <= 19

	subject to cons2:
       u + x(1) >= 0
       u + x(2) >= -1
       u + x(3) >= -2
       u + x(4) >= -3
       u + x(5) >= -4
       u + x(6) >= -5
       u + x(7) >= -6
       u + x(8) >= -7
       u + x(9) >= -8
       u + x(10) >= -9
       u + x(11) >= -10
       u + x(12) >= -11
       u + x(13) >= -12
       u + x(14) >= -13
       u + x(15) >= -14
       u + x(16) >= -15
       u + x(17) >= -16
       u + x(18) >= -17
       u + x(19) >= -18
       u + x(20) >= -19

問題を解きます。

[sol,fval,exitflag,output] = solve(prob)
Solving problem using quadprog.

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

<stopping criteria details>

sol = struct with fields:
    u: 4.1786
    x: [20×1 double]


fval = 653.3036

exitflag = 
    OptimalSolution


output = struct with fields:
            message: 'Minimum found that satisfies the constraints.↵↵Optimization completed because the objective function is non-decreasing in ↵feasible directions, to within the value of the optimality tolerance,↵and constraints are satisfied to within the value of the constraint tolerance.↵↵<stopping criteria details>↵↵Optimization completed: The relative dual feasibility, 1.421085e-16,↵is less than options.OptimalityTolerance = 1.000000e-08, the complementarity measure,↵1.536130e-10, is less than options.OptimalityTolerance, and the relative maximum constraint↵violation, 0.000000e+00, is less than options.ConstraintTolerance = 1.000000e-08.'
          algorithm: 'interior-point-convex'
      firstorderopt: 6.1445e-09
    constrviolation: 0
         iterations: 8
       linearsolver: 'sparse'
       cgiterations: []
             solver: 'quadprog'

解での変数 x を表します。

disp(sol.x)
    1.5000
    3.0000
    4.5000
    6.0000
    7.5000
    9.0000
   10.1786
   11.1786
   12.1786
   13.1786
   14.1786
   15.1786
   16.1786
   17.1786
   18.1786
   19.1786
   20.1786
   21.1786
   22.1786
   23.1786

解での制約関数値を表示します。

evaluate(cons1,sol)
ans = 1×20

   -2.6786   -2.1786   -1.6786   -1.1786   -0.6786   -0.1786   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000


evaluate(cons2,sol)
ans = 1×20

    5.6786    8.1786   10.6786   13.1786   15.6786   18.1786   20.3571   22.3571   24.3571   26.3571   28.3571   30.3571   32.3571   34.3571   36.3571   38.3571   40.3571   42.3571   44.3571   46.3571

制約 cons1 は、値が非正であることを指定します。cons1 の最初の 6 つの値は負であり、残りの値は精度を表示するために 0 となります。そのため、解では、最初の 6 つの値を除くすべての値について制約 cons1 がアクティブとなります。これに対し、制約 cons2 は値が非負であることを指定します。cons2 のすべての値は厳密に正であり、この制約が解においてアクティブでないことを意味します。

補助関数

次のコードは、forloopfcn 補助関数を作成します。

function [cons1,cons2] = forloopfcn(x,u)
N = length(x);
cons1 = zeros(1,N);
cons2 = zeros(1,N);
for i = 1:N
    cons1(i) = x(i) - u - i + 1;
    cons2(i) = x(i) + u + i - 1;
end
end

参考

トピック