Niniejsza praca proponuje ramy obliczeń fizycznych łączące wspólny Temporal Exponential Random Graph Model z ulepszonym algorytmem Motter-Lai w celu oceny ewolucji strukturalnej i ilościowego określenia progu odporności ekosystemów gospodarki cyfrowej w warunkach kaskadowych awarii.
Artykuł badawczy
Ewolucja strukturalna i odporność ekosystemów gospodarki cyfrowej: podejście oparte na połączeniu czasowego wykładniczego modelu losowych grafów oraz ram teoretycznych Mottera-Laia
W tym artykule
Podsumowanie
Streszczenie
Istniejące metody analizy statycznej pomijają ewolucję strukturalną złożonych topologii sieci oraz kaskadowe awarie wywołane lokalną redystrybucją obciążenia, co prowadzi do błędów w ocenie odporności ekosystemów gospodarki cyfrowej. Aby dokładnie określić próg odporności systemu, w niniejszej pracy zaproponowano ramy obliczeń fizycznych łączące wspólny czasowy wykładniczy model losowych grafów (Temporal Exponential Random Graph Model – TERGM) z ulepszonym algorytmem Mottera-Laia (TERGM-ML). Ramy te wykorzystują estymację maksymalnej wiarygodności metodą Monte Carlo dla łańcuchów Markowa (MCMC-MLE) do modelowania endogennych efektów strukturalnych i rekonstrukcji trajektorii ewolucji czasowej rzeczywistej topologii sieci, co pozwala przezwyciężyć ograniczenia modeli statycznych. Następnie, w oparciu o centralność węzłów i ich nieliniową wydajność fizyczną, w przypadku celowego ataku uruchamiana jest reguła redystrybucji ruchu zależna od pozostałej wydajności sąsiadów, co umożliwia śledzenie całego procesu dezintegracji systemu spowodowanego propagacją lokalnego przeciążenia. Symulacje porównawcze wielu modeli wykazują, że po wprowadzeniu podwójnego mechanizmu ewolucji czasowej i dynamicznej realokacji, krytyczny próg usuwania węzłów wyzwalający załamanie globalnej wydajności transmisji w scenariuszu celowego ataku opartego na centralności pośrednictwa wynosi 12,41% ± 0,63%, co jest wartością znacznie wyższą niż w statycznym modelu bazowym sieci bezskalowej (7,85% ± 0,42%, p < 0,001).
Wprowadzenie
Wraz z głęboką integracją globalnej technologii informacyjnej, ekosystem gospodarki cyfrowej ewoluował stopniowo w złożony system sieciowy o charakterze transgranicznym, splecionym i wysoce współzależnym1,2. Badania nad ewolucją strukturalną i granicami odporności tego systemu mają ogromne znaczenie strategiczne dla zapewnienia stabilnego funkcjonowania makroekonomii oraz bezpieczeństwa przemysłu cyfrowego. Jednakże, w obliczu przekształceń globalnego krajobrazu gospodarczego i częstych asymetrycznych wstrząsów zewnętrznych, widoczna stała się podatność topologiczna sieci wykazana przez ekosystem gospodarki cyfrowej3,4. Istniejące ekonometryczne metody makroekonomiczne i statyczne analizy statystyczne często traktują wewnętrzne relacje systemu jako liniowe kombinacje zmiennych, co uniemożliwia ujawnienie dynamicznych praw przejścia fazowego w odpowiedzi na ryzyka ekstremalne z perspektywy mikrotopologicznego kaskadowania5,6,7.
Aby zaradzić tym ograniczeniom, opracowano fizyczny model obliczeniowy z połączonym czasowym wykładniczym modelem losowych grafów (TERGM)8,9 oraz ulepszonym algorytmem Mottera-Laia, co pozwala rozwiązać problemy techniczne związane z dyskretnością czasową, założeniami dotyczącymi podziału obciążenia oraz rozdzielnością architektur w istniejących badaniach10. Model ten wypełnia luki matematyczne między interakcjami na poziomie mikro a awarią całego systemu na poziomie makro poprzez pomiar granicy odporności ekosystemu gospodarki cyfrowej poddanego planowanym atakom asymetrycznym. Zintegrowany model obliczeniowy opracowany w niniejszej pracy nie tylko usprawnia mechanizm wnioskowania matematycznego w zakresie odporności ewolucji sieci złożonych, ale także dostarcza wysoce powtarzalnych podstaw matematycznych do zapobiegania kryzysom związanym z globalnymi zakłóceniami sieci w erze cyfrowej.
Protokół
Protokół składa się z czterech kolejnych etapów obliczeniowych, które przekształcają empiryczne dane panelowe w ilościowy próg odporności dla ekosystemu gospodarki cyfrowej.
Ewolucja topologii czasowej za pomocą TERGM
Pierwszy etap polega na wykonaniu ewolucji topologii przy użyciu modelu czasowego eksponencjalnego grafu losowego z estymacją maksymalnego prawdopodobieństwa opartą na metodzie Monte Carlo w łańcuchu Markowa. Empiryczny zbiór danych panelowych ICT-DE500, reprezentujący podmioty gospodarki cyfrowej w latach 2018–okres obserwacji 2023, jest importowany do środowiska systemowego, a numery indeksów jednostek są wyrównywane w różnych krokach czasowych w celu utworzenia sekwencji grafów skierowanych pionowo. Wielkość kroku aktualizacji parametrów jest ustalona na 0,01, a pierwsze 10 000 iteracji MCMC są odrzucane jako okres wstępnego rozgrzewania (burn-in), aby osiągnąć rozkład stacjonarny. Zestaw danych ICT-DE500 został zbudowany na podstawie globalnych rekordów inwestycji korporacyjnych i fuzji z Crunchbase obejmujących lata 2018–2023, przy czym jednostki zostały przefiltrowane tak, aby uwzględniać wyłącznie podmioty działające w sektorach Technologii Informacyjno-Komunikacyjnych oraz Gospodarki Cyfrowej. W tej skonstruowanej sieci węzły reprezentują pojedyncze podmioty gospodarki cyfrowej (tj. przedsiębiorstwa i instytucje inwestycyjne), a krawędzie skierowane reprezentują nieważone binarne przepływy kapitałowe poprzez inwestycje lub fuzje i przejęcia (M&A)&A) zdarzenia. Te interakcje finansowe stanowią logiczne ścieżki propagacji modelowanego obciążenia, ponieważ zależności kapitałowe i przepływy finansowe tworzą bezpośrednie kanały transmisji ryzyka; trudności finansowe w jednym węźle wymuszają ponowne rozłożenie płynności i deprecjację aktywów, które bezpośrednio przenoszą się na topologicznie połączone podmioty. 500 podstawowych podmiotów wybrano na podstawie najwyższego poziomu centralności w sieci oraz najaktywniejszych rejestrów interakcji w okresie obserwacji. Pokrojono krawędzie roczne dla każdego z sześciu lat. W celu zapewnienia ścisłej synchronizacji czasowej numerów indeksów podmiotów oraz ujednolicenia wymiarów macierzy (N = 500) dla estymacji TERGM, odosobnione węzły o stopniu zero w danym roku'płaty s zostały zachowane jako tymczasowo nieaktywne jednostki, a nie usunięte strukturalnie. Zbieżność TERGM ocenia się poprzez monitorowanie trajektorii parametrów MCMC-MLE dla wszystkich endogenicznych współczynników strukturalnych, w tym gęstości krawędzi, wzajemności oraz geometrycznie ważonej liczby wspólnych partnerstw krawędziowych. Łańcuch uznaje się za zbieżny, gdy wszystkie trajektorie parametrów wykazują stabilne oscylacje wokół odpowiednich średnich wartości bez kierunkowego dryfu po przekroczeniu progu 10 000 kroków okresu wypalania (burn-in). Po osiągnięciu zbieżności parametrów system wykonuje 10 000 kolejnych iteracji próbkowania Gibbsa w celu modelowania efektów endogenicznych, takich jak tendencje do tworzenia centrów gwiazdowych, generując ciągłe, czasowo zsynchronizowane topologie reprezentujące sieć'makroskopowej ewolucji strukturalnej
Aby sformalizować proces generowania, specyfikacja matematyczna modelu TERGM modeluje prawdopodobieństwo warunkowe zaobserwowania docelowej topologii sieci Gt w kroku czasowym makroskopowym t, przy założonej poprzedniej sieci Gt-1, ponieważ
Tutaj, θ jest wektorem parametrów podstawowych kontrolującym ewolucję strukturalną, h(Gt, Gt-1) to sieć's wektorem statystyk wystarczalnych ilościowo opisujących wcześniej opisane składniki strukturalne endogenne (tj. gęstość krawędzi, wzajemność oraz geometrycznie ważone wspólne powiązania krawędziowe), oraz c(θ, Gt-1to funkcja podziału zapewniająca normalizację prawdopodobieństwa. Dla kolejnych symulacji nieliniowych kaskadowych uszkodzeń, końcowa stabilna realizacja sieci pochodzi z wygenerowanego poprzedniego ciągu G1:T jest wyodrębniane, aby służyć jako początkowy substrat topologiczny. Zasadniczo, ponieważ ewolucja strukturalna makroekonomiczna przebiega w znacznie dłuższej skali czasowej (makrokroki czasowe, t) niż chwilowe lokalne awarie kaskadowe, topologia sieci nie ewoluuje dalej za pomocą mechanizmów TERGM podczas symulacji kaskady. Zamiast tego zmiany topologiczne zachodzące w trakcie szybkich kroków czasowych mikrokaskady (τsą wywoływane wyłącznie celowym usuwaniem węzłów oraz powstającymi w wyniku przeciążenia awariami wtórnymi.
Kalibracja pojemności fizycznej i inicjalizacja obciążenia
Drugi etap wykonuje kalibrację pojemności fizycznej dla wszystkich węzłów w sekwencji macierzy topologii sieci wyjściowej. Dla każdego węzła wyodrębniane są całkowity stopień oraz skierowalna centralność pośrednictwa, z małą stałą o wartości 10-8 wprowadzony w obliczeniach centralności pośrednictwa, aby uniknąć dzielenia przez zero spowodowanego dyskretnością lokalnej sieci. Początkowe obciążenie usługi Li(0) jest mapowane na wszystkie węzły sieci za pomocą nieliniowego równania potęgowego
gdzie ki to znormalizowany całkowity stopień, Bi to znormalizowana skierowana centralność pośrednictwa, λ to współczynnik wagi równowagi (ustalony na 0,5, aby zapewnić równą wagę), oraz β to zakres ograniczeń indeksu alokacji obciążenia wynosi od 1,0 do 1,5. Limit pojemności fizycznej przenoszenia Ci dla każdego węzła jest ustalana poprzez zastosowanie hiperparametru tolerancji pojemności na poziomie systemu α (od 0,1 do 0,5) w celu utworzenia fizycznej granicy odporności na wstrząsy
Dolne ograniczenie α = 0,1 oznacza minimalny scenariusz nadmiarowości, w którym węzły posiadają jedynie 10% zapasowej pojemności powyżej swojego podstawowego obciążenia, podczas gdy górna granica α = 0,5 odpowiada konfiguracji o wysokiej nadmiarowości z 50% zapasowym zapotrzebowaniem. Wartości pośrednie α = 0,2, 0,3 i 0,4 są również stosowane w analizie wrażliwości dwuzmiennej w celu skonstruowania pełnej ortogonalnej przestrzeni parametrów z wykładnikiem heterogeniczności obciążenia βIndeks alokacji obciążenia β mieści się w zakresie od 1,0 do 1,5, gdzie β = 1,0 daje liniowy rozkład obciążenia i β = 1,5 generuje silnie spolaryzowane skupienie obciążenia w kierunku węzłów o wysokiej centralności. Współczynnik wagi równowagi λ jest ustalony na 0,5, aby zapewnić równy wkład stopnia i centralności pośrednictwa w obliczeniach początkowego obciążenia. Ustawienia parametrów rdzenia dla symulacji ewolucji szeregów czasowych i kaskadowych uszkodzeń podsumowano w Tabela 1.
Nieliniowa dynamika kaskadowa w przypadku zamierzonych ataków
Trzeci etap implementuje nieliniową dynamikę kaskadową w warunkach zamierzonego ataku. Symulacja uruchamia zamierzony atak poprzez blokadę i wymuszone usunięcie zestawu węzłów centralnych w kolejności malejącej centralności pośrednictwa, co narusza początkową ochronę topologii cyfrowego ekosystemu w celu przetestowania dynamicznego obciążenia w mikrokrokach kaskadowych. Silnik przeładunku obciążenia jest uruchamiany w celu ponownego kierowania przepływu nadmiarowego, ograniczonego przez rzeczywistą, chwilową pojemność fizyczną węzłów sąsiednich, z członem ujścia o wartości 10-8 wprowadzony w celu symulacji przepełnienia aktywów cyfrowych, gdy drogi komercyjne są całkowicie zablokowane. Węzeł uznawany jest za uszkodzony, gdy jego obciążenie chwilowe przekracza jego pojemność fizyczną, a weryfikacja tego przeciążenia wykonywana jest równolegle we wszystkich aktywnych węzłach w celu aktualizacji binarnej funkcji stanu przeżycia. Zamierzony atak skierowany jest na węzły w ściśle malejącej kolejności centralności pośrednictwa, przy czym każdy krok ataku usuwa dokładnie jeden węzeł ze zbioru przetrwających węzłów aktywnych. Waga przeprofilowania obciążenia 
przypisane z węzła, który uległ awarii i ∈ Fτ do sąsiada, który przeżył j ∈ Aτ w kroku czasowym mikro τ jest obliczane jako
gdzie 
reprezentuje pozostałą pojemność fizyczną sąsiada j, Gij jest wskaźnikiem topologicznej przyległości, Aτ jest aktywnym zbiorem węzłów przetrwałych, a 
= 10-8 zapobiega dzieleniu przez zero. Na podstawie tych wag, skala obciążenia chwilowego 
liczba przetrwałych węzłów jest synchronicznie zmieniana
Następnie weryfikuje się aktualizację stanu awarii wtórnej za pomocą binarnej funkcji przeżycia 
:
Węzeł jest deklarowany jako nieprawidłowy (
) gdy jego obciążenie chwilowe przekracza pojemność, aktualizacja zestawu uszkodzeń Fτ+1. Kaskada osiąga stan ustalony, gdy Fτ+n = ∅, co wskazuje, że w bieżącej mikrokroku czasowym nie doszło do dodatkowych uszkodzeń węzłów i wszystkie przetrwające węzły działają w granicach swoich możliwości.
Cykl kaskadowy trwa, aż żadne dodatkowe węzły nie ulegną awarii, co oznacza osiągnięcie wtórnego stanu ustalonego, w którym potencjał kaskady został całkowicie rozproszony.
Ocena odporności systemu i identyfikacja progów
Czwarty etap ocenia odporność systemu poprzez monitorowanie makroskopowego wskaźnika rozpadu grafu cyfrowego ekosystemu. W sposób ciągły wyodrębnia się względną skalę największej spójnej składowej powstałej z węzłów, które przetrwały, w celu sporządzenia krzywej rozpadu przejścia fazowego w funkcji frakcji usuniętych węzłów. Globalna efektywność transmisji E(τ) jest obliczane w celu ilościowego określenia spójności przetrwałej topologii
gdzie N to początkowa całkowita liczba węzłów (stała), Aτ to zbiór przetrwających aktywnych węzłów, oraz 
to najkrótsza skierowana odległość geodezyjna od węzła i do j na bieżącym etapie. Próg przejścia krytycznego jest następnie identyfikowany poprzez monitorowanie zmian pierwszej pochodnej tej funkcji wydajności względem współczynnika usunięcia. Próg ten wyznacza się jako punkt, w którym pierwsza pochodna osiąga wartość minimalną, co wskazuje na najbardziej stromy spadek wydajności transmisji. Próg odporności krytycznej oblicza się poprzez różniczkowanie numeryczne globalnej wydajności transmisji E(τ) względem współczynnika usunięcia węzłów f przy użyciu drugiego centralnego schematu różnicowego rzędu. Tor pochodnej pierwszego rzędu dE/df jest wygładzany przy pomocy średniej ruchomej z pięciu kolejnych punktów danych, aby zmniejszyć szum próbkowania Monte Carlo, zachowując jednocześnie położenie najbardziej stromego spadku. Próg krytyczny fc jest wybierany jako współczynnik usunięcia, przy którym wygładzona pierwsza pochodna osiąga swoją wartość globalnego minimum, odpowiadającą punktowi maksymalnej szybkości spadku sprawności transmisji. Kryterium to jest stosowane konsekwentnie we wszystkich scenariuszach symulacji oraz modelach podstawowych. Podana wartość progowa wynosi 12,41% ± 0,63% reprezentuje średnią i odchylenie standardowe obliczone z 100 niezależnych symulacji Monte Carlo z różnymi ziarnami losowymi, zapewniając statystyczną odporność lokalizacji przejścia fazowego.
Konfiguracje symulacji i implementacje podstawowe
Aby zapewnić pełną powtarzalność symulacji, przed każdą iteracją Monte Carlo przypisywano kolejne wartości ziarna losowego (liczby całkowite od 1 do 100). Ewolucję topologiczną i modelowanie statystyczne wykonano w języku R z wykorzystaniem pakietu tergm, natomiast nieliniowe symulacje kaskadowe zaimplementowano w języku Python przy użyciu biblioteki NetworkX. Dodatkowo, dla porównawczego modelu głębokiego uczenia wykorzystano model GCN-Attack zaimplementowany w PyTorch Geometric. Zbudowano go z wykorzystaniem standardowej dwuwarstwowej architektury grafowej sieci konwolucyjnej (ukryta wymiarowość 64) i trenowano przy użyciu optymalizatora Adam z szybkością uczenia 0,01 przez 200 epok, aby zapewnić rygorystyczny i spójny kontrolny przebieg eksperymentu podczas oceny modeli referencyjnych.
Wyniki
Ogólna logika działania i przepływ danych proponowanego ramach obliczeń fizycznych są zilustrowane w Rycina 1. W miarę działania struktury framework są rejestrowane mikroskopowe cechy termiczne lokalnego przeładunku i przełączania obciążenia oraz nieliniowa ewolucja rozkładów stopnia węzłów (omówione w Rycina 2 i Rycina 3, a dynamiczne szczegóły opisano poniżej. Kolejne sekcje bezpośrednio odwzorowują wyniki symulacji na etapy protokołu.
Ewolucja topologii czasowej za pomocą TERGM
Rycina 4 wizualnie przedstawia dekonstrukcję przestrzennej topologii i struktury społecznościowej rdzeniowej sieci ICT-DE500, podkreślając rozmieszczenie węzłów o wysokiej centralności pośrednictwa, które były celem symulowanych zamierzonych ataków. Test dobroci dopasowania potwierdza, że wygenerowana topologia sieci skutecznie modeluje ewolucję czasową rzeczywistych ekosystemów, skutecznie unikając eksplozji gradientów lub pułapek optymalności lokalnej po okresie rozgrzewki składającym się z 10 000 kroków. Rycina 5 przedstawia trajektorie diagnostyczne zbieżności parametrów MCMC-MLE oraz dobrze dopasowania wyrażone odległością geodezyjną. Rycina 5A pokazuje, że trzy podstawowe parametry reprezentujące gęstość krawędzi θ₁, wzajemność θ₂oraz geometrycznie ważona wspólna partycypacja krawędziowa θ₃ wszystkie przerywają swój duży, kierunkowy dryf po przekroczeniu progu burn-in wynoszącego 10 000 kroków, przy czym oczekiwane średnie zbliżają się do poziomej linii bazowej i stabilizują się w jej pobliżu. Rycina 5B demonstruje, że empiryczne obserwacje najkrótszych odległości geodezyjnych zdecydowanie mieszczą się w granicach rozkładu uzyskanego dla 1000 niezależnych realizacji sieci. Realizacje te wyodrębniono z 10 000 kolejnych iteracji próbkowania Gibbsa, stosując odstęp rzadzenia równy 10, aby zminimalizować autokorelację, co potwierdza wiarygodność podstawy generowania topologii. Szczegółowe oszacowania parametrów MCMC-MLE, błędy standardowe oraz istotność statystyczna efektów endogenicznych strukturalnych w poszczególnych latach obserwacji są przedstawione w Tabela 2.
Ewolucja czasowa makroskopowej struktury topologicznej jest ilościowo określana w Rysunek 6Gęstość sieci wzrastała systematycznie z poziomu 0,015 do 0,035 w latach 2018–2023, podczas gdy średni współczynnik grupowania wzrósł z 0,22 do 0,37, co wskazuje na istotne zjawisko rozbieżności między gęstością a grupowaniem. Najbardziej stromy wzrost gęstości odnotowano w latach 2020–2021, kiedy wzrosła ona z 0,021 do 0,029, natomiast współczynnik grupowania osiągnął lokalne maksimum wynoszące około 0,31 w 2020 roku, by następnie spaść do około 0,29 mimo szybkiego wzrostu gęstości w 2021 roku. Ta rozbieżność ujawnia mechanizm adaptacyjnej ewolucji w warunkach fluktuacji cyklu makroekonomicznego, w którym grupowanie związane z unikaniem ryzyka w 2020 roku sprzyja lokalnemu zagęszczaniu, podczas gdy masowe nowe połączenia transgraniczne w 2021 roku tymczasowo osłabiają strukturę ścisłych społeczności.
Kalibracja pojemności fizycznej i inicjalizacja obciążenia
Analiza czułości dwuzmiennowa w Rysunek 7 analizuje łączny wpływ nadmiarowości pojemności fizycznej i polaryzacji obciążenia na trajektorię przejścia fazowego największego spójnego komponentu. W obrębie dziewięciu ortogonalnych kombinacji tolerancji pojemności α i heterogeniczność obciążenia β, zestaw paneli pokazuje, że zwiększenie α i zmniejszanie się β obie opóźniają kolaps sieci. W scenariuszu obciążenia spolaryzowanego z β = 1,5 i minimalna redundancja α = 0,1 w Rycina 7A, próg krytycznego załamania wynosi około fc = 0,08. Podnoszenie α do 0,5 w Rycina 7C przesuwa punkt przegięcia w prawo do fc ≈ 0,23. W warunkach obciążenia zrównoważonego z β = 1,0 i α = 0,1 w Rycina 7G, próg pozostaje odporny przy fc ≈ 0,18 oraz przy optymalnym połączeniu α = 0,5 i β = 1,0 w Rycina 7I, próg znacząco wydłuża się do fc ≈ 0,38. Wyniki te wykazują, że równoważenie obciążenia przynosi większy korzyści marginalne w zakresie odporności niż samoumnażanie pojemności.
Nieliniowa dynamika kaskadowa w warunkach zamierzonych ataków
Jak przedstawiono we wprowadzeniu do ramowego przeglądu, mikroskopowe cechy termiczne lokalnej przebudowy obciążenia po początkowym awaryjnym kaskadzie są pokazane w Rysunek 2, a nieliniowa ewolucja rozkładu stopnia węzłów w trzech typowych mikrokrokach czasowych została przedstawiona w Rysunek 3.
Ocena odporności systemu i identyfikacja progów
Próg przejścia krytycznego dla globalnej efektywności transmisji znajduje się na poziomie 12,41% ± 0,63% usunięcie węzłów podczas ataku ukierunkowanego. W kontekście sieci o 500 węzłach, ten udział odpowiada ukierunkowanemu usunięciu około 62 węzłów centralnych. Ten próg oznacza punkt załamania efektywności (tj. początek najbardziej stromego spadku efektywności transmisji), a nie całkowite rozłączenie topologiczne. Rycina 8 przedstawia trójwymiarową powierzchnię ewolucji globalnej wydajności E(τ) nad współczynnikiem usuwania i szczytowym obciążeniem sieci Rycina 8A, a dwuwymiarowy przekrój z numerycznym różniczkowaniem w Rycina 8B. Gdy współczynnik usunięcia f jest poniżej 0,10, E(τ) pozostaje powyżej 0,8, a pierwsza pochodna oscyluje w obszarze płytkim. Minimum trajektorii pierwszej pochodnej wskazuje próg przejścia krytycznego, przy czym Rycina 8B wyświetlanie pojedynczego przekroju poprzecznego przy częstotliwości fc = 12,0%, co jest wysoce zgodne ze średnią statystyczną z 100 niezależnych symulacji Monte Carlo.
Konfiguracje symulacji i implementacje podstawowe
Proponowany model znacząco przewyższa modele statyczne i modele głębokiego uczenia w scenariuszach ataków ukierunkowanych. Jednak w warunkach awarii losowych statyczny model BA-ML wykazuje wyższy próg przeżycia (49,12%) w porównaniu z modelem TERGM-ML (46,28%). Należy podkreślić, że porównanie z modelem statycznym BA-ML stanowi odrębną topologiczną belkę pomiarową, a nie ściśle kontrolowane usunięcie, ponieważ model BarabáMechanizm generatywny si-Alberta istotnie różni się od ramy ERGM. Rycina 9 wyświetla wykres „deszczowo-chmurowy” szczytowych prędkości propagacji kaskady w czterech architekturach modeli. Model bazowy Static BA-ML wykazuje medianę szczytowej prędkości wynoszącą około 49,7 węzła na krok, przy ekstremalnych przypadkach zbliżających się do 140. Modele SNA-Cascading i GCN-Attack mają mediany odpowiednio około 35,6 i 23,9. Model TERGM-ML wykazuje najsilniejszą zbieżność z medianą 13,2 węzła na krok, niemal całkowicie eliminując ekstremalne kolapsy przekraczające 40. Tabela 3 podsumowuje krytyczne progi i istotność statystyczną dla wszystkich modeli. Te porównania jasno pokazują, że choć ramy TERGM-ML wykazują lepszą odporność strukturalną na ukierunkowane asymetryczne szoki, zaobserwowane różnice w wydajności wynikają z łącznego wpływu różnych podstawowych topologii, ewolucji czasowej oraz przeprojektowania z uwzględnieniem pojemności, a nie wyłącznie z pojedynczych wyłączeń mechanizmów.
DOSTĘPNOŚĆ DANYCH:
Nieprzetworzone dane wykorzystane w tym badaniu pochodzą z globalnej bazy danych inwestycji korporacyjnych i fuzji Crunchbase, dostępnej publicznie na platformie Kaggle pod adresem https://www.kaggle.com/datasets/justinas/startup-investments. Przetworzony podzbiór ICT-DE500, składający się z 500 podmiotów z rocznymi macierzami krawędzi za okres 2018 roku–Dane z 2023 roku oraz dane atrybutów węzłów, w tym stopień i centralność pośrednictwa, wraz z kodami szacowania modelu TERGM i diagnostyki zbieżności, kodem symulacji kaskadowego uszkodzenia z ulepszonym algorytmem Mottera-Lai oraz pełnymi specyfikacjami zależności, zostały umieszczone w publicznym repozytorium GitHub pod adresem https://github.com/moonmoon1189/digital-economy-resilience-complex-networks.
Rycina 1: Ewolucja topologii czasowej i nieliniowa kaskadowa struktura obliczeń fizycznych. Ten rysunek przedstawia ogólną logikę wykonywania i przepływ danych, w tym etapy ewolucji topologii, kalibracji pojemności fizycznej, nieliniowego kaskadowania oraz oceny odporności, służące identyfikacji progu przejścia krytycznego. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.
Rycina 2: Diagram cieplny mikroskopowej ewolucji nieliniowego przepływu obciążenia i lokalnego przepełnienia kaskadowego. Rysunek przedstawia dynamiczne cechy termiczne lokalnej redistribucji obciążenia przepływu po początkowym uszkodzeniu kaskadowym, od mikrokroku czasowego 0 do kroku 5. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.
Rysunek 3Wykres rozrzutu nieliniowej ewolucji rozkładu stopnia węzłów podczas awarii kaskadowej. Na rysunku przedstawiono trajektorię ewolucji rozkładu stopnia węzłów systemu w trzech charakterystycznych mikrokrokach czasowych (0, 3, 6). Kliknij tutaj, aby wyświetlić większą wersję tego rysunku.
Rycina 4: Topologia przestrzenna, struktura społeczności oraz rozkład celów ataków ukierunkowanych w sieci rdzeniowej ICT-DE500. Rysunek przedstawia wizualną dekompozycję wysoce nieliniowej topologii makroskopowej oraz mikroskopowych atrybutów węzłów sieci rdzennej, z zaznaczeniem centrów typu gwiazda i podatnych na uszkodzenia źródeł. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.
Rycina 5Test zgodności trajektorii diagnostycznych parametrów metodą Monte Carlo łańcucha Markowa oraz odległości geodezyjnej. (A) Ten panel pokazuje trajektorię diagnostyczną estymacji parametrów MCMC-MLE w kolejnych iteracjach, podczas gdy panel (Bwyświetla test dobroci dopasowania najkrótszej odległości geodezyjnej. Kliknij tutaj, aby wyświetlić większą wersję tego rysunku.
Rycina 6: Wspólny rozkład parametrów makrotopologicznych cech w ewolucji szeregów czasowych. Na rysunku przedstawiono zmieniające się trendy makrotopolologicznych parametrów, a mianowicie gęstości sieci oraz średniego współczynnika grupowania, dla ekosystemu gospodarki cyfrowej w latach 2018–2023. Zakresy cieniowane wokół linii trendu odpowiadają 95% przedziałom ufności uzyskanym na podstawie 100 niezależnych symulacji Monte Carlo. Kliknij tutaj, aby wyświetlić większą wersję tego rysunku.
Rysunek 7Wykres tablicowy rozpadu przejścia fazowego dla dwuzmiennej wrażliwości na tolerancję pojemności i heterogeniczność obciążenia. (A–ITe panele przedstawiają trajektorie przejść fazowych dla różnych ortogonalnych kombinacji tolerancji pojemności i heterogeniczności obciążenia. Kliknij tutaj, aby wyświetlić większą wersję tego rysunku.
Rysunek 8: Wspólny profil globalnej wydajności transmisji w trzech wymiarach oraz progowy przekrój krytycznej zmiany Panel (A) konstruuje trójwymiarową ewolucję przestrzenną globalnej sprawności transmisyjnej, a panel (Bwyodrębnia lokalizację dwuwymiarowego progowego przejścia krytycznego (tj. punktu załamania sprawności) w przekroju poprzecznym, stosując różniczkowanie numeryczne. Kliknij tutaj, aby wyświetlić większą wersję tego rysunku.
Rysunek 9: Niejednorodny rozkład prędkości szczytowego rozprzestrzeniania się kaskady w wykresie chmurowym Rysunek kompleksowo przedstawia heterogeniczny rozkład gęstości prawdopodobieństwa prędkości szczytowych propagacji kaskady czterech modeli podczas wybuchów katastrof wtórnych. Kliknij tutaj, aby wyświetlić większą wersję tego rysunku.
| Parametr | Zmienna & Zakres | Ograniczenie & Uzasadnienie |
| Skala sieci | N=500 | Stała wielkość wyodrębnionego podzbioru empirycznej sieci ICT-DE500. |
| Waga centralności | λ=0.5 | Ukorzeniona linia bazowa zapewniająca równą wagę stopniowi i pośrednictwu. |
| Wskaźnik alokacji obciążenia | β∈[1.0,1.5] | Nieliniowy wykładnik kontrolujący początkową heterogeniczność obciążenia biznesowego. |
| Tolerancja pojemnościowa | α∈[1.0,1.5] | Stosunek nadmiarowości na poziomie systemu określający górne ograniczenie pojemności fizycznej. |
| Okres wypalenia MCMC | >10 000 kroków | Obowiązkowe początkowe iteracje odrzucone w celu osiągnięcia rozkładu stacjonarnego. |
| Iteracje próbkowania Gibbsa | >10 000 kroków | Kolejne kroki próbkowania umożliwiające uzyskanie zgodnych czasowo topologii sieci. |
Tabela 1: Podstawowe ustawienia parametrów dla symulacji ewolucji szeregów czasowych i kaskadowych uszkodzeń fizycznych. W tabeli zdefiniowano podstawowe parametry, w tym skalę sieci, wagę centralności, indeks alokacji obciążenia, tolerancję pojemności oraz iteracje metody Monte Carlo łańcucha Markowa (MCMC).
| Rok obserwacji | Gęstość krawędzi (θ1) Szacowanie | Gęstość krawędzi (θ1) Błąd std. | Współwystępowanie (θ2) Szacowanie | Współwystępowanie (θ2) Błąd std. | GWESP (θ3) Oszacowanie | GWESP (θ3) Błąd std. | Znaczenie |
| 2018 | -6.35 | 0.12 | 1.95 | 0.08 | 1.25 | 0.05 | *** |
| 2019 | -6.15 | 0.11 | 1.98 | 0.09 | 1.35 | 0.06 | *** |
| 2020 | -6.05 | 0.14 | 2.15 | 0.1 | 1.55 | 0.07 | *** |
| 2021 | -5.85 | 0.13 | 1.9 | 0.09 | 1.3 | 0.06 | *** |
| 2022 | -5.75 | 0.12 | 2.05 | 0.08 | 1.42 | 0.05 | *** |
| 2023 | -5.65 | 0.11 | 2.1 | 0.07 | 1.48 | 0.05 | *** |
Tabela 2: Szacunki parametrów modelu TERGM dla endogenicznych efektów strukturalnych w poszczególnych latach obserwacji (2018–2023). W tabeli przedstawiono oszacowania parametrów, błędy standardowe oraz istotność statystyczną uzyskane za pomocą estymacji największej wiarygodności w oparciu o łańcuch Markowa Monte Carlo (MCMC-MLE) dla gęstości krawędzi, wzajemności oraz geometrycznie ważonej liczby wspólnych partnerstw krawędziowych w sześciu rocznych okresach obserwacji. ***p < 0,001. Błędy standardowe podano obok szacunków parametrów.
| Architektura modelu | Mechanizm ewolucji czasowej | Mechanizm dynamicznego przeznaczania | Próg krytyczny (atak ukierunkowany) | Próg krytyczny (awaria losowa) | Istotność statystyczna (wartość p) |
| TERGM-ML | Tak | Tak | 12.41% ± 0.63% | 46.28% ± 1.75% | Odniesienie – linia bazowa |
| Statyczny BA-ML | Nie | Tak | 7.85% ± 0.42% | 49.12% ± 1.88% | p < 0.001 *** |
| SNA-Kaskadowanie | Tak | Nie | 8.93% ± 0.55% | 37.54% ± 1.42% | p = 0,003 ** |
| GCN-Attack (SOTA Baseline) | Niejawny | Niejawny | 10.76% ± 0.81% | 43.15% ± 2.05% | p = 0,021 * |
Tabela 3: Porównanie granicy odporności mechanizmu rdzeniowego ablacji i architektury wielomodelowej. Tabela zawiera szczegółowe informacje o krytycznych progach oraz wynikach testów statystycznych odporności systemu dla grafów szeregów czasowych i ulepszonego modelu Motter-Lai (TERGM-ML) oraz trzech modeli bazowych w scenariuszach zamierzonych ataków i losowych awarii. Wartości podano jako średnią ± odchylenie standardowe na podstawie 100 niezależnych symulacji Monte Carlo. Atak ukierunkowany oznacza sekwencyjne usuwanie węzłów na podstawie malejącej centralności pośrednictwa. Istotność statystyczna ocenia różnicę w progu ataku ukierunkowanego między odpowiednim modelem bazowym a proponowanym podejściem przy użyciu niezależnego testu t dla dwóch prób (*p < 0.05, **p < 0.01, ***p < 0.001).
Dyskusja
Zaproponowane kaskadowe ramy obliczeń fizycznych łączące wspólne grafy szeregów czasowych i ulepszony model Mottera-Laia (TERGM-ML) skutecznie niwelują ograniczenia typu „czarna skrzynka” tradycyjnych modeli opartych wyłącznie na danych w prognozowaniu odporności. Ramy te opierają się na fundamentalnych wykładniczych modelach losowych grafów wprowadzonych przez Wassermana i Pattisona11 oraz na ramowej strukturze ataków kaskadowych opracowanej pierwotnie przez Mottera i Laia12, rozszerzając oba podejścia o dynamikę temporalną i lokalne ograniczenia przepustowości. Paradygmat ten ściśle wiąże rzeczywistą endogenną ewolucję topologiczną z limitami obciążenia mikroelementów poprzez wprowadzenie selektywnej logiki przepływu opartej na lokalnych ograniczeniach granic przepustowości fizycznej. Mechanizm ewolucji temporalnej jest zgodny ze specyfikacjami TERGM dla dynamicznego modelowania sieci13,14, a strategia alokacji przepustowości jest spójna z zasadami projektowania redundancji sieci w celu łagodzenia awarii kaskadowych15,16.
Kluczowym krokiem w protokole jest optymalny mechanizm przekierowania podstawowego przepływu biznesowego w oparciu o dostępną wydajność sąsiadów, zastępujący nierealistyczne założenie o „średnim rozkładzie” w tradycyjnym modelu Mottera-Laia. Założenie o jednolitym redystrybuowaniu w standardowym modelu Mottera-Laia było krytykowane w niedawnych badaniach nad odpornością infrastruktury za pominięcie heterogenicznych ograniczeń wydajności węzłów6,10. Niniejsze wyniki wskazują, że podstawowa sieciowa struktura społeczności w pętli zamkniętej wykazuje wyraźny fizyczny efekt tłumienia szczytów przeciążeń, skutecznie hamując propagację kaskadową i znacznie opóźniając przejście fazowe prowadzące do dezintegracji globalnej wydajności transmisji. Model TERGM-ML posiada najwyższy próg krytyczny dla celowych ataków, osiągający 12,41% ± 0,63%, co odzwierciedla zdolność tłumienia endogennej architektury sieci i łagodzi globalne ryzyko lawinowe wywołane pojedynczym punktem przeciążenia. Podniesienie progu krytycznego z 7,85% do 12,41% wynika z dwóch synergicznych mechanizmów. Mechanizm ewolucji czasowej generuje struktury społeczności w pętli zamkniętej oraz więzi wzajemne, których brakuje w statycznych sieciach bezskalowych. Społeczności te ograniczają przestrzennie propagację przeciążeń, wymuszając na nadmiarowym obciążeniu przejście przez wiele ścieżek wewnątrz społeczności przed dotarciem do odległych obszarów, przy czym każdy krok przejścia rozprasza część obciążenia przejściowego poprzez absorpcję przez sąsiednie węzły. Mechanizm dynamicznej redystrybucji kieruje nadmiarowe obciążenie wyłącznie do sąsiadów z dodatnią pozostałą wydajnością ΔCj(τ) > 0, unikając jednolitego rozkładu, który w standardowym modelu Mottera-Laia szybko wyczerpuje lokalną nadmiarowość. Społeczności w pętli zamkniętej zapewniają strukturę topologiczną, która sprawia, że routing świadomy wydajności jest skuteczny, podczas gdy routing świadomy wydajności zapobiega przedwczesnemu nasyceniu łączy wewnątrz społeczności. To sprzężenie wyjaśnia, dlaczego zintegrowany model przewyższa statyczną linię bazową o ponad 4 punkty procentowe w progu krytycznym. Wartość tego progu jest zgodna z teoretycznymi przewidywaniami dla sieci bezskalowych poddawanych ukierunkowanym atakom17 oraz z zachowaniami przejścia fazowego perkolacji obserwowanymi w systemach złożonych18.
Pomimo tych osiągnięć, metoda ta posiada pewne ograniczenia. Ze względu na istniejące granice obserwacyjne, obecne ekstrapolacje w dużym stopniu opierają się na scentralizowanych, kompletnych przekrojach globalnej topologii, a ich dyskretne okna czasowe próbkowania nie pozwalają na dokładne uchwycenie mikrozmiennych w czasie zaburzeń impedancji wywołanych przez wysokoczęstotliwościowe, nagłe oscylacje środowiska zewnętrznego. Ograniczenia te odzwierciedlają wyzwania zidentyfikowane w niedawnych przeglądach wskaźników odporności dla systemów cyber-fizycznych oraz modelowania kaskadowych awarii w warunkach dynamicznych19,20. Przyszłe badania i zastosowania mogą zostać rozszerzone o architektury zdecentralizowane, koncentrując się na badaniu adaptacyjnych mechanizmów dynamicznej kompensacji odporności w oparciu o rozproszoną współpracę wieloagentową w warunkach gier z niepełną informacją. Abstrakcja grafu jednowarstwowego oraz globalne przypisanie parametrów stanowią kluczowe ograniczenia obecnego modelu. Badania nad sieciami wielowarstwowymi wykazały, że współzależności między warstwami interakcji mogą wzmacniać lub tłumić propagację kaskady w sposób, którego modele jednowarstwowe nie są w stanie uchwycić. Globalne przypisanie tolerancji pojemności α oraz wykładnika alokacji obciążenia β pomija specyficzną dla poszczególnych podmiotów heterogeniczność marginesów pojemności i czułości obciążenia. W przyszłych pracach warto zbadać trzy kierunki rozszerzenia: zastąpienie topologii jednowarstwowej reprezentacją wielowarstwową, która rozróżnia przepływy kapitałowe, licencjonowanie technologii i świadczenie usług jako oddzielne warstwy z zależnościami międzywarstwowymi; kalibrację parametrów pojemności i obciążenia specyficznych dla podmiotów na podstawie operacyjnych danych na poziomie przedsiębiorstw; oraz przejście od scentralizowanych przekrojów topologii do zdecentralizowanych architektur wieloagentowych, w których węzły podejmują adaptacyjne decyzje o redystrybucji w oparciu o lokalnie obserwowalne sygnały. Niedawne badania nad sieciami wielowarstwowymi wykazały, że współzależności między odrębnymi warstwami interakcji mogą wzmacniać lub tłumić propagację kaskady w sposób, którego modele jednowarstwowe nie są w stanie uchwycić.
Oświadczenia
Autorzy oświadczają, że nie istnieją żadne konflikty interesów. W procesie tworzenia, generowania lub modyfikowania jakichkolwiek elementów graficznych nie wykorzystano narzędzi generatywnej sztucznej inteligencji (AI).
Wkład autorów:
F.Y. i Y.Z. opracowali koncepcję i zaprojektowali badanie. F.Y. przeprowadził symulacje obliczeniowe, przeanalizował dane i przygotował pierwotny szkic manuskryptu. Y.Z. nadzorował badania, zapewnił wsparcie teoretyczne oraz krytycznie zrewidował manuskrypt pod kątem istotnej zawartości intelektualnej. Wszyscy autorzy zapoznali się z ostateczną wersją manuskryptu i ją zatwierdzili.
Podziękowania
Autorzy nie otrzymali wsparcia od żadnej organizacji w związku z przedłożoną pracą.
Materiały
| Nazwa | Firma | Numer katalogowy | Komentarze |
|---|---|---|---|
| AMD EPYC 7742 CPU | Advanced Micro Devices | 7742 | Wydajny procesor do przeszukiwania struktur grafów i przeliczania najkrótszych ścieżek. |
| Crunchbase Database | Kaggle | startup-investments | Globalne rekordy inwestycji korporacyjnych oraz sieci fuzji i przejęć (M&A) wykorzystane jako baza sieci globalnej. |
| CUDA 11.6 | NVIDIA | version 11.6 | Platforma akceleracji sprzętowej wykorzystywana do operacji tensorowych w modelu bazowym GCN. |
| NetworkX 2.8 | NetworkX Developers | version 2.8 | Biblioteka do analizy złożonych sieci używana do ekstrakcji parametrów grafu i wyszukiwania ścieżek. |
| NumPy | NumPy Developers | N/A | Framework jądra matematycznego zapewniający deterministyczną logikę i eliminujący dryft numeryczny. |
| NVIDIA RTX 3090 GPU | NVIDIA | RTX 3090 | Procesor graficzny wdrożony w celu akceleracji obliczeń tensorowych w modelu bazowym głębokiego uczenia. |
| Python 3.9 | Python Software Foundation | version 3.9 | Podstawowe środowisko wykonawcze, w którym skompilowano i uruchomiono główny framework. |
| PyTorch 1.12 | Meta AI | version 1.12 | Biblioteka głębokiego uczenia używana do obliczeń grafowych i propagacji w przód w modelu bazowym. |
| R | / | version 4.2.2 | / |
| statnet package | The statnet Project | N/A | Zaawansowany pakiet rozszerzeń statystycznych używany do wieloetapowego TERGM MCMC-MLE dla sieci dynamicznych. |
| tergm package | / | version 4.2.0 | / |
| Ubuntu 22.04.1 LTS | Canonical | 22.04.1 LTS | Konfiguracja systemu operacyjnego serwera hostująca wielowątkową macierz obliczeniową. |
Bibliografia
- Rong K. Research agenda for the digital economy. J Digit Econ. 2022;1(1):20–31.
- Fan R, et al. Network dynamics of inter-firm innovation in China’s digital economy: a two-layer network perspective. Technol Anal Strateg Manag. 2025:1–20.
- Feng Y, Huang M. The geographical analysis of global economic uncertainty: resource distribution, geopolitical risks, and systemic vulnerability. Geogr Res Bull. 2025;4:570–573.
- Zhang H, Liu H, Chen R. Multilayer innovation network resilience: a framework for digital economy vulnerability assessment. iScience. 2026;29(1):114295.
- Zang T, et al. Current status and perspective of vulnerability assessment of cyber-physical power systems based on complex network theory. Energies. 2023;16(18):6509.
- He S, et al. Cascading failure in cyber-physical systems: a review on failure modeling and vulnerability analysis. IEEE Trans Cybern. 2024;54(12):7936–7954.
- Dong G, Sun Z, Sun N, Wang F. Understanding percolation phase transition behaviors in complex networks from the macro and meso-micro perspectives. Europhys Lett. 2022;139(6):61001.
- Shi X, Huang X, Liu H. Research on the structural features and influence mechanism of the low-carbon technology cooperation network based on temporal exponential random graph model. Sustainability. 2022;14(19):12341.
- Yao X, Du Y, Pu Y, Wang B. Structural evolution and its determinants of domestic value-added network of digital service exports based on temporal exponential random graph model. Emerg Mark Finance Trade. 2024;60(14):3387–3401.
- Lu Z, Qiu W. Resilience analysis of seaport-dry-port network in container transport: multi-stage load redistribution dynamics following cascade failure. Systems. 2025;13(4):299.
- Wasserman S, Pattison P. Logit models and logistic regressions for social networks: I. An introduction to Markov graphs and p*. Psychometrika. 1996;61(3):401–425.
- Motter AE, Lai YC. Cascade-based attacks on complex networks. Phys Rev E. 2002;66(6):065102.
- Fritz C, Mehrl M, Thurner PW, Kauermann G. Exponential random graph models for dynamic signed networks: an application to international relations. Polit Anal. 2025;33(3):211–230.
- Li Y, Pu Y. Pattern evolution and dynamic formation mechanism of global scrap copper trade network: based on temporal exponential random graph model. Ecol Econ. 2025;236:108664.
- Liu J, Liu X, Liu P. Capacity allocation strategy against cascading failure of complex network. J Syst Eng Electron. 2024;35(6):1507–1515.
- Motter AE. Cascade control and defense in complex networks. Phys Rev Lett. 2004;93(9):098701.
- Albert R, Jeong H, Barabási AL. Error and attack tolerance of complex networks. Nature. 2000;406(6794):378–382.
- Artime O, et al. Robustness and resilience of complex networks. Nat Rev Phys. 2024;6(2):114–131.
- Li ZS, Wu G, Cassandro R, Wang H. A review of resilience metrics and modeling methods for cyber-physical power systems. IEEE Trans Reliab. 2024;73(1):59–66.
- Ma C, et al. A review of supply chain resilience: a network modeling perspective. Appl Sci. 2025;15(1):265.