ペトリネットと状態方程式


※記事自体を差し替えました
おまけにペトリネットジェネレーターをつけました

ペトリネットのモデルは以下の形の式で表されます

$x(k+1) = x(k) + Bu(k) $


この式の形は離散値系の状態方程式

$x(k+1) = Ax(k) + Bu(k) $
$y(k) = Cx(k) + Du(k) $

において、$A$が単位行列の特別な場合と解釈できます。
出力$y$はおいておきます

つまり、スケジュリーング問題において制御理論が適用可能となります。

この辺はペトリネットを用いた以下の論文においても言及されています。
Philipp Wenzelburger ∗ Frank Allg¨ower  
"A Petri Net Modeling Framework for the Control of Flexible Manufacturing Systems"  

https://www.ist.uni-stuttgart.de/institute/team/pdf/PW/MIM19_PN_Modeling.pdf  



具体例として以下の論文に簡単な実装例が乗っています。
Kenji Sawada*, Seiichi Shin*, Kenji Kumagai** and Hisato Yoneda**
"Optimal Scheduling of Automatic Guided Vehicle System via State Space Realization"

https://www.fujipress.jp/ijat/au/ijate000700050571/

私が実装した際はソルバとしてGLPKを採用しました。理由としては、論文中の1つ目のアルゴリズムにおいて
初期値のまま使用した際に解が得られた唯一のソルバでした。(他:Matlab,Cbc,CPLEXで試しました.Mathematicaは未実施)


おまけとして、最初に挙げた論文に記載されているalgorithm1～3を用いたペトリネット生成器の実装例を貼ります。
なお、論文内に一部抜粋で記載されている状態方程式(3)式は間違っています。

B行列の4列目、入力$[m_1,m_2,\tau_3,\tau_5,o_3,S]$に対する列ですが、
$\begin{matrix}[0 &-1& 0& 0& 0& 0&  1& 0& 0& 0 ]\^T \end{matrix}$ となっており、状態変数の定義と照らし合わせると
$[m_1 \tau_1,o_3,B] から  [m_2 \tau_3,o_3,B]$の遷移を示しているため、3列目と取り違えていると思われます。
入力の記法によると,$m_1\tau_3からm_2\tau_5$への遷移のはずなので。

以下実装例　記載のままほぼベタ書きですのですぐ理解いただけると思います。英語やアルゴリズムの誤りがあったらご指摘いただけると非常に嬉しいです。　
```julia
# Julia v1.7.2
using LinearAlgebra

task = 5
order = 3
machine = 2
kp = ones(Int8,machine,task) # kp matrix is production steps at each pairs of machine and task.
kp[1,3] = 3 

tm = [1 2 3; # able to do tasks in machine. row 1 means machine No.1, row 2 machine No.2 ...
      3 4 5];
to = [1 2 3; # task of order
      1 3 4;
      1 3 5];

p = String[] # places
tʳ = String[] # transitions
x = Vector{Int8}[] # state value

# Algorithm 1. Create Places 53 places
for m = 1:machine
    push!(p,"p(𝑚$(m),0,0,𝐼)")
end
for o = 1:order
    push!(p,"p(0,0,𝑜$(o),𝑆)")
    for t in to[o,:]
        push!(p,"p(0,τ$(t),𝑜$(o),𝐶)")
        push!(p,"p(0,τ$(t),𝑜$(o),𝑁)")
        for m = 1:machine
            if t ∈ tm[m,:]
               for production = 1:(kp[m,t])
                    push!(p,"p(𝑚$(m),τ$(t),𝑜$(o),𝑃₍$(production)₎)")
               end
               push!(p,"p(𝑚$(m),τ$(t),𝑜$(o),𝐵)")
            end
        end
    end
end

# 依存関係説明関数 t を始めるまでに終わっている必要のあるtask
function necessaryToStartTheTask(t)
    if t == 2 || t == 3 || t == 4
        return [1]
    elseif t == 5
        return [3]
    else
        return []
    end
end

# Algorithm 2. Determine indirect dependencies
function fn(t, t̂)
    # t1 being necesarry to start any others
    # t3 being necesarry to start t5
    for t̄ in necessaryToStartTheTask(t̂)
        if t̄ == t
            return 1
        else
            return fn(t,t̄)
        end
    end
    return 0
end

# Algorithm 3. Create Transitions and Arcs
A = Matrix{Int8}(I,length(p),length(p))
B = Matrix{Int8}(undef,length(p),0)
isEmptySet = true
for o = 1:3#1:order
    for t in to[o,:]
        for m = 1:machine
            if t ∈ tm[m,:]
                # if task t does not require any prior task
                if t == 1
                    # transition
                    push!(tʳ,"t(0,𝑚$(m),0,τ$(t), 𝑜$(o),𝑆)")
                    
                    # input arcs
                    column1 = findfirst(p.=="p(0,τ$(t),𝑜$(o),𝑁)")
                    column2 = findfirst(p.=="p(𝑚$(m),0,0,𝐼)")
                    column3 = findfirst(p.=="p(0,0,𝑜$(o),𝑆)")
                    # output arc
                    row1 = findfirst(p.=="p(𝑚$(m),τ$(t),𝑜$(o),𝑃₍1₎)")
                    #A[row1,column1] = 1
                    #A[row1,column2] = 1
                    #A[row1,column3] = 1
                    bVec = zeros(Int8,length(p))
                    bVec[column1] = -1
                    bVec[column2] = -1
                    bVec[column3] = -1
                    bVec[row1] = 1
                    global B  
                    B = hcat(B,bVec)
                end
        
                # production transition
                for production = 1:(kp[m,t]-1)
                    # transition
                    push!(tʳ,"t(𝑚$(m),𝑚$(m),τ$(t),τ$(t), 𝑜$(o),𝑃₍$(production)₎)")
                    # input arcs
                    column1 = findfirst(p.=="p(𝑚$(m),τ$(t),𝑜$(o),𝑃₍$(production)₎)")
                    A[column1,column1] = 0 #対角項の削除
                    # output arc
                    row1 = findfirst(p.=="p(𝑚$(m),τ$(t),𝑜$(o),𝑃₍$(production+1)₎)")
                    A[row1,column1] = 1
                end

                # Finishing transition
                push!(tʳ,"t(𝑚$(m),𝑚$(m),τ$(t),τ$(t), 𝑜$(o),𝐹)")
                column1 = findfirst(p.=="p(𝑚$(m),τ$(t),𝑜$(o),𝑃₍$(kp[m,t])₎)")
                row1 = findfirst(p.=="p(𝑚$(m),τ$(t),𝑜$(o),𝐵)")
                row2 = findfirst(p.=="p(𝑚$(m),0,0,𝐼)")
                row3 = findfirst(p.=="p(0,τ$(t),𝑜$(o),𝐶)")
                A[column1,column1] = 0 #対角項の削除
                A[row1,column1] = 1
                A[row2,column1] = 1
                A[row3,column1] = 1

                # ever task t ∈ to
                for t̂ in to[o,:]
                    for m̂ = 1:machine
                        if t̂ ∈ tm[m̂,:]
                        # 以下論文の条件の通りfn(t,tᵟ) == fn(tᵟ,t̂) == 1を満たすtᵟ ∈ toの集合が零集合であること
                            if t != t̂ && fn(t̂,t) == 0
                                global isEmptySet = true
                                for tᵟ in to[o,:]
                                    if fn(t,tᵟ) == 1 && fn(tᵟ,t̂) == 1
                                        global  isEmptySet = false
                                    end
                                end
                                if isEmptySet
                                    # if task t does not require any prior task
                                    # transition
                                    push!(tʳ,"t(𝑚$(m),𝑚$(m̂),τ$(t),τ$(t̂), 𝑜$(o),𝑆)")
                                    # input arcs
                                    column1 = findfirst(p.=="p(0,τ$(t̂),𝑜$(o),𝑁)")
                                    column2 = findfirst(p.=="p(𝑚$(m̂),0,0,𝐼)")
                                    column3 = findfirst(p.=="p(𝑚$(m),τ$(t),𝑜$(o),𝐵)")
                                    # output arc
                                    row1 = findfirst(p.=="p(𝑚$(m̂),τ$(t̂),𝑜$(o),𝑃₍1₎)")
                                    bVec = zeros(Int8,length(p))
                                    bVec[column1] = -1
                                    bVec[column2] = -1
                                    bVec[column3] = -1
                                    bVec[row1] = 1
                                    for tΔ in necessaryToStartTheTask(t̂)
                                        # transition
                                        column1 = findfirst(p.=="p(0,τ$(tΔ),𝑜$(o),𝐶)")
                                        row1 = findfirst(p.=="p(0,τ$(tΔ),𝑜$(o),𝐶)")
                                        
                                        #A[row1,column1] = 1
                                        bVec[column1] = -1
                                        bVec[row1] = 1
                                    end
                                    global B  
                                    B = hcat(B,bVec) # ここのtransitionは1つなのでbVecは１つ

                                end
                            end
                        end
                    end
                end
            end
        end
    end
end
# 変数名はUnicodeで遊んでいただけです。
🦄 = unique(tʳ); # tʳに重複がある場合があったため重複を削除
# 以下、論文内(3),(4)式と比べるための確認部分
cutA1 = findfirst(p.=="p(𝑚1,τ1,𝑜3,𝑃₍1₎)")
cutA2 = findfirst(p.=="p(𝑚2,τ5,𝑜3,𝐵)")
🦄start = 🦄[findall(x ->occursin("𝑆",x),🦄)]
m = 1
t = 1
o = 3
cutB1 = findfirst(🦄.=="t(0,𝑚$(m),0,τ$(t), 𝑜$(o),𝑆)")
cutB11 = findfirst(🦄start.=="t(0,𝑚$(m),0,τ$(t), 𝑜$(o),𝑆)")
m = 2
m̂ = 2
t = 3
t̂ = 5
o = 3
cutB2 = findfirst(🦄.=="t(𝑚$(m),𝑚$(m̂),τ$(t),τ$(t̂), 𝑜$(o),𝑆)")
cutB22 = findfirst(🦄start.=="t(𝑚$(m),𝑚$(m̂),τ$(t),τ$(t̂), 𝑜$(o),𝑆)")
𝑥ᴵᴰ= p[cutA1:cutA2,:]
sampleA = A[cutA1:cutA2,cutA1:cutA2]
sampleA = sampleA[[1,2,5,6,7,8,9,10,13,14],[1,2,5,6,7,8,9,10,13,14]]
𝑢ᴵᴰ =  🦄[cutB1:cutB2]
sampleB = B[cutA1:cutA2,cutB11:cutB22]
sampleB = sampleB[[1,2,5,6,7,8,9,10,13,14],:]
```