Если SPiCT — это гипотеза о состоянии запаса, а HCR — рецепт, то MSE — испытательный стенд, на котором рецепт готовят в турбулентной кухне реальности. Это «аэродинамическая труба» для управления: мы создаём правдоподобный мир (операционную модель), пропускаем через него шум наблюдений, оценку, задержки внедрения, ошибки исполнения — и смотрим, как ведут себя разные стратегии в тысячах альтернативных будущих. Не чтобы найти «идеально оптимальную» — её нет, — а чтобы выбрать наименее хрупкую.
Из чего состоит MSE.
- Операционная модель (OM): «истинная» динамика запаса и флота, с вариабельностью рекрутинга, автокорреляцией, возможным дрейфом уловистости, изменчивостью среды.
- Наблюдательная модель: как рождаются данные (ошибки индексов, задержки отчётности, недоучёт).
- Оценка (например, SPiCT): периодическая, с неизбежными смещениями и сходимостями.
- Исполнение: реализация с ошибками (implementation error), реакция флота, соблюдение ограничений, лаги решений.
Что считаем успехом.
- Экономика/социальная стабильность: средний вылов, межгодовая изменчивость ОДУ, доля «закрытых» лет, длительность восстановления.
- Компромисс: фронт Парето вместо одного числа; подход «минимакс сожаления» — минимум проигрыша в плохих сценариях.
Чем MSE спасает от иллюзий.
- Мы перестаём путать модель с реальностью: оценки систематически ошибаются, и это встроено в эксперименты.
- Широкие сценарии важнее глубины одной «идеальной» модели: узкие допущения делают стратегию хрупкой.
- Поведение людей имеет значение: соблюдение квот, перераспределение усилия, экономические стимулы — часть петли обратной связи.
Что делаем в этой главе.
- Формулируем цели и метрики: явно, измеримо, с порогами риска.
- Задаём наблюдательный и оценочный контур (в т. ч. SPiCT), календарь пересмотров, лаги и ошибки исполнения.
- Сравниваем HCR‑кандидаты, настраиваем их параметры под целевые риски, считаем фронт компромиссов.
- Проводим стресс‑тесты: «плохие годы подряд», резкий технологический скачок, смещение индекса, неполное соблюдение.
- Готовим «панель управления» для решений: таблицы производительности, графики риска/вылова/стабильности, понятные без кода.
Границы и дисциплина. MSE не отвечает на вопрос «что истина?» — он отвечает «что работает чаще и не ломается в плохую погоду». Он требует прозрачного кода, фиксированных зерён случайности и версионирования сценариев. И, как хорошие инженерные мосты, лучшие MSE‑решения реже всего «самые умные» — они просто переживают больше штормов.
# =============================================================== # СКРИПТ 3: MANAGEMENT STRATEGY EVALUATION (MSE)- ОЦЕНКА СТРАТЕГИИ УПРАВЛЕНИЯ # Сравнение трех ключевых стратегий управления # Курс: Оценка водных биоресурсов при недостатке данных в R # Автор: Баканёв С.В. # Дата создания: 28.08.2025 # =============================================================== # ======================= ВВЕДЕНИЕ ============================= # В этом скрипте сравниваются три основные стратегии управления: # 1. Fish at Fmsy - промысел на уровне оптимальной смертности # 2. MSY hockey-stick rule - адаптивное правило с защитой запаса # 3. ICES advice rule - предосторожный подход ICES # # MSE позволяет тестировать эти стратегии с учетом всех # источников неопределенности в системе промысел-запас # ------------------- 1. ПОДГОТОВКА СРЕДЫ -------------------- ## 1.1 Очистка рабочей среды и настройка rm (list = ls ())set.seed (123 ) # Для воспроизводимости результатов ## 1.2 Загрузка необходимых библиотек library (spict) # Для работы с моделью SPiCT
Загрузка требуемого пакета: TMB
Welcome to spict_v1.3.8@107a32
library (tidyverse) # Обработка данных
-- Attaching core tidyverse packages ------------------------ tidyverse 2.0.0 --
v dplyr 1.1.4 v readr 2.1.5
v forcats 1.0.0 v stringr 1.5.1
v ggplot2 3.5.2 v tibble 3.2.1
v lubridate 1.9.4 v tidyr 1.3.1
v purrr 1.0.4
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag()
i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library (ggplot2) # Визуализация library (viridis) # Цветовые схемы
Загрузка требуемого пакета: viridisLite
library (patchwork) # Компоновка графиков library (gridExtra) # Дополнительные возможности компоновки
Присоединяю пакет: 'gridExtra'
Следующий объект скрыт от 'package:dplyr':
combine
## 1.3 Установка рабочей директории и загрузка модели setwd ("C:/SPICT" )# Загружаем подогнанную модель из первого скрипта if (file.exists ("spict_model_fit.rds" )) {<- readRDS ("spict_model_fit.rds" )cat (" \n ✓ Модель SPiCT успешно загружена \n " )else {stop ("Файл модели не найден. Запустите первый скрипт для подгонки модели." )
<U+2713> Модель SPiCT успешно загружена
cat (" \n " , strrep ("=" , 60 ), " \n " )
============================================================
cat (" MANAGEMENT STRATEGY EVALUATION (MSE) \n " )
MANAGEMENT STRATEGY EVALUATION (MSE)
cat (" Сравнение трех стратегий управления промыслом \n " )
Сравнение трех стратегий управления промыслом
cat (strrep ("=" , 60 ), " \n " )
============================================================
# ------------------- 2. ИЗВЛЕЧЕНИЕ ПАРАМЕТРОВ ИЗ МОДЕЛИ -------------------- cat (" \n ========== ПАРАМЕТРЫ ОПЕРАЦИОННОЙ МОДЕЛИ ========== \n " )
========== ПАРАМЕТРЫ ОПЕРАЦИОННОЙ МОДЕЛИ ==========
## 2.1 Извлечение оценок параметров популяционной динамики # Эти параметры представляют "истинное" состояние в операционной модели <- get.par ("logr" , fit, exp = TRUE )[1 ] # Внутренний темп роста <- get.par ("logK" , fit, exp = TRUE )[1 ] # Ёмкость среды <- get.par ("logBmsy" , fit, exp = TRUE )[1 ] # Биомасса MSY <- get.par ("logFmsy" , fit, exp = TRUE )[1 ] # Промысловая смертность MSY <- get.par ("logB" , fit, exp = TRUE )[1 ] # Текущая биомасса <- get.par ("logF" , fit, exp = TRUE )[1 ] # Текущая пром. смертность # Параметр улавливаемости для индекса CPUE <- get.par ("logq" , fit, exp = TRUE )[1 ]## 2.2 Вывод ключевых параметров cat (" \n Извлеченные параметры модели: \n " )
Извлеченные параметры модели:
cat (sprintf ("r (темп роста): %.3f год⁻¹ \n " , r_true))
r (темп роста): 0.289 год<U+207B><U+00B9>
cat (sprintf ("K (ёмкость средыь): %.1f тыс. т \n " , K_true))
K (ёмкость средыь): 153.8 тыс. т
cat (sprintf ("Bmsy: %.1f тыс. т \n " , Bmsy_true))cat (sprintf ("Fmsy: %.3f год⁻¹ \n " , Fmsy_true))
Fmsy: 0.145 год<U+207B><U+00B9>
cat (sprintf ("Текущая биомасса: %.1f тыс. т \n " , B_current))
Текущая биомасса: 103.0 тыс. т
cat (sprintf ("Текущее B/Bmsy: %.2f \n " , B_current/ Bmsy_true))cat (sprintf ("Текущее F/Fmsy: %.2f \n " , F_current/ Fmsy_true))# ------------------- 3. НАСТРОЙКИ СИМУЛЯЦИИ -------------------- cat (" \n ========== НАСТРОЙКИ MSE ========== \n " )
========== НАСТРОЙКИ MSE ==========
## 3.1 Основные параметры симуляции <- 500 # Количество симуляций (реализаций) <- 100 # Период прогнозирования (лет) <- 2 # Интервал между оценками запаса (лет) cat (sprintf ("Количество симуляций: %d \n " , n_sim))
Количество симуляций: 500
cat (sprintf ("Горизонт прогнозирования: %d лет \n " , n_years))
Горизонт прогнозирования: 100 лет
cat (sprintf ("Частота оценки запаса: каждые %d года \n " , assessment_interval))
Частота оценки запаса: каждые 2 года
## 3.2 Параметры неопределенности # Все источники неопределенности в системе # Процессная ошибка - естественная изменчивость в динамике популяции <- 0.15 # Ошибка наблюдения - неточность в индексах биомассы <- 0.20 # Ошибка оценки - неточность в оценке состояния запаса <- 0.15 # Ошибка реализации - разница между рекомендованным и фактическим выловом <- 0.10 cat (" \n Параметры неопределенности: \n " )
Параметры неопределенности:
cat (sprintf ("CV процессной ошибки: %.0f%% \n " , process_error_cv * 100 ))
CV процессной ошибки: 15%
cat (sprintf ("CV ошибки наблюдения: %.0f%% \n " , observation_error_cv * 100 ))
CV ошибки наблюдения: 20%
cat (sprintf ("CV смещения в оценках: %.0f%% \n " , assessment_bias_cv * 100 ))
CV смещения в оценках: 15%
cat (sprintf ("CV ошибки реализации: %.0f%% \n " , implementation_error_cv * 100 ))
CV ошибки реализации: 10%
# ------------------- 4. ОПЕРАЦИОННАЯ МОДЕЛЬ -------------------- cat (" \n ========== ОПРЕДЕЛЕНИЕ ОПЕРАЦИОННОЙ МОДЕЛИ ========== \n " )
========== ОПРЕДЕЛЕНИЕ ОПЕРАЦИОННОЙ МОДЕЛИ ==========
## 4.1 Функция истинной динамики популяции # Модель Шефера с процессной стохастичностью <- function (B_t, F_t, r, K, process_cv = 0.15 ) {# Расчет прибавочной продукции (модель Шефера) <- r * B_t * (1 - B_t/ K)# Расчет вылова <- F_t * B_t# Обновление биомассы с учетом процессной ошибки # Используем лог-нормальное распределение для мультипликативной ошибки <- rlnorm (1 , meanlog = - process_cv^ 2 / 2 , sdlog = process_cv)<- (B_t + surplus_production - catch) * process_error# Ограничения для реалистичности <- max (B_next, 0.001 * K) # Минимум 0.1% от K <- min (B_next, 1.5 * K) # Максимум 150% от K return (list (B_next = B_next,catch_realized = catch,surplus = surplus_production,process_multiplier = process_error## 4.2 Функция генерации наблюдаемого индекса # Имитирует процесс сбора данных с ошибками <- function (B_true, q, obs_cv = 0.20 ) {# Истинный индекс пропорционален биомассе <- q * B_true# Добавление ошибки наблюдения <- rlnorm (1 , meanlog = - obs_cv^ 2 / 2 , sdlog = obs_cv)<- true_index * obs_errorreturn (observed_index)## 4.3 Функция оценки состояния запаса # Упрощенная процедура оценки (в реальности здесь бы запускалась полная модель) <- function (index_history, catch_history, assessment_cv = 0.15 ) {# Оценка текущей биомассы по последним наблюдениям индекса # Используем среднее за последние 3 года для сглаживания <- tail (index_history[! is.na (index_history)], 3 )<- mean (recent_indices, na.rm = TRUE )# Оценка биомассы с учетом смещения <- rlnorm (1 , meanlog = - assessment_cv^ 2 / 2 , sdlog = assessment_cv)<- (mean_recent_index / true_q) * assessment_bias# Оценка текущей F из последнего вылова <- tail (catch_history[catch_history > 0 ], 1 )<- ifelse (length (recent_catch) > 0 && B_estimated > 0 ,/ B_estimated,0.1 ) # Значение по умолчанию # Оценки референсных точек также содержат неопределенность <- true_Bmsy * rlnorm (1 , 0 , assessment_cv/ 2 )<- true_Fmsy * rlnorm (1 , 0 , assessment_cv/ 2 )return (list (B = B_estimated,F = F_estimated,Bmsy = Bmsy_estimated,Fmsy = Fmsy_estimated,B_Bmsy = B_estimated / Bmsy_estimated,F_Fmsy = F_estimated / Fmsy_estimated# ------------------- 5. ПРАВИЛА УПРАВЛЕНИЯ (HCR) -------------------- cat (" \n ========== ОПРЕДЕЛЕНИЕ ПРАВИЛ УПРАВЛЕНИЯ ========== \n " )
========== ОПРЕДЕЛЕНИЕ ПРАВИЛ УПРАВЛЕНИЯ ==========
## 5.1 Стратегия 1: Fish at Fmsy # Простейшее правило - всегда промысел на уровне Fmsy <- function (assessment, previous_TAC = NULL ) {# Рекомендация по F <- assessment$ Fmsy# Расчет TAC <- F_advice * assessment$ B# Обеспечиваем неотрицательность <- max (TAC, 0 )return (list (F_advice = F_advice,TAC = TAC,rule_name = "Fish at Fmsy" ,status = ifelse (assessment$ B_Bmsy < 0.5 , "Риск истощения" , "Стандартный промысел" )cat ("✓ Стратегия 1: Fish at Fmsy \n " )
<U+2713> Стратегия 1: Fish at Fmsy
cat (" - Постоянный промысел на уровне Fmsy \n " )
- Постоянный промысел на уровне Fmsy
cat (" - Не учитывает состояние запаса \n " )
- Не учитывает состояние запаса
cat (" - Простое в применении правило \n " )
- Простое в применении правило
## 5.2 Стратегия 2: MSY Hockey-stick Rule # Адаптивное правило с защитой при низкой биомассе <- function (assessment, previous_TAC = NULL ,Blim_fraction = 0.5 ) {# Определение предельной биомассы <- assessment$ Bmsy * Blim_fraction# Применение правила hockey-stick if (assessment$ B <= Blim) {# Линейное снижение F при B < Blim <- assessment$ B / Blim<- assessment$ Fmsy * F_multiplier<- "Снижение промысла" else {# Полный промысел при B >= Blim <- assessment$ Fmsy<- "Полный промысел" # Расчет TAC <- F_advice * assessment$ B<- max (TAC, 0 )# Полное закрытие при критически низкой биомассе if (assessment$ B_Bmsy < 0.2 ) {<- 0 <- 0 <- "Промысел закрыт" return (list (F_advice = F_advice,TAC = TAC,rule_name = "MSY Hockey-stick" ,status = status,Blim = Blimcat (" \n ✓ Стратегия 2: MSY Hockey-stick Rule \n " )
<U+2713> Стратегия 2: MSY Hockey-stick Rule
cat (" - Адаптивное управление в зависимости от состояния \n " )
- Адаптивное управление в зависимости от состояния
cat (" - Линейное снижение F при B < 0.5*Bmsy \n " )
- Линейное снижение F при B < 0.5*Bmsy
cat (" - Полное закрытие при B < 0.2*Bmsy \n " )
- Полное закрытие при B < 0.2*Bmsy
## 5.3 Стратегия 3: ICES Advice Rule # Предосторожный подход с ограничениями межгодовых изменений <- function (assessment, previous_TAC = NULL ,Bpa_multiplier = 1.4 ,Fpa_multiplier = 0.85 ,max_TAC_change = 0.20 ) {# Расчет предосторожных ориентиров управления <- assessment$ Bmsy / Bpa_multiplier # Предосторожная биомасса <- Bpa / 1.4 # Предельная биомасса <- assessment$ Fmsy * Fpa_multiplier # Предосторожная F # Определение F по правилу ICES if (assessment$ B < Blim) {# Критическое состояние - закрытие промысла <- 0 <- "Критическое - промысел закрыт" else if (assessment$ B >= Blim && assessment$ B < Bpa) {# Восстановление - пропорциональное снижение F <- (assessment$ B - Blim) / (Bpa - Blim)<- Fpa * F_multiplier<- "Восстановление запаса" else {# Нормальное состояние - предосторожный промысел <- Fpa<- "Устойчивый промысел" # Расчет TAC <- F_advice * assessment$ B<- max (TAC, 0 )# Ограничение межгодовых изменений TAC (стабильность для промышленности) if (! is.null (previous_TAC) && previous_TAC > 0 && TAC > 0 ) {<- previous_TAC * (1 + max_TAC_change)<- previous_TAC * (1 - max_TAC_change)<- min (max (TAC, max_decrease), max_increase)if (TAC != TAC_constrained) {<- paste (status, "(TAC ограничен)" )<- TAC_constrainedreturn (list (F_advice = F_advice,TAC = TAC,rule_name = "ICES Advice Rule" ,status = status,Bpa = Bpa,Blim = Blim,Fpa = Fpacat (" \n ✓ Стратегия 3: ICES Advice Rule \n " )
<U+2713> Стратегия 3: ICES Advice Rule
cat (" - Предосторожный подход (Fpa = 0.85*Fmsy) \n " )
- Предосторожный подход (Fpa = 0.85*Fmsy)
cat (" - Двухуровневая система пороговых значений \n " )
- Двухуровневая система пороговых значений
cat (" - Ограничение межгодовых изменений TAC (±20%) \n " )
- Ограничение межгодовых изменений TAC (±20%)
# ------------------- 6. ФУНКЦИЯ MSE ДЛЯ ОДНОЙ СИМУЛЯЦИИ -------------------- ## 6.1 Основная функция симуляции <- function (sim_id, HCR_function, params, settings) {# Инициализация с небольшой вариацией начальных условий <- params$ B_initial * rlnorm (1 , 0 , 0.05 )# Вариация в "истинных" параметрах (представляет неопределенность в природе) <- params$ r * rlnorm (1 , 0 , 0.05 )<- params$ K * rlnorm (1 , 0 , 0.05 )# Массивы для хранения результатов <- settings$ n_years<- data.frame (sim_id = sim_id,year = 1 : n_years,B_true = numeric (n_years),B_estimated = numeric (n_years),F_advice = numeric (n_years),F_realized = numeric (n_years),catch = numeric (n_years),TAC = numeric (n_years),index_obs = numeric (n_years),B_Bmsy_true = numeric (n_years),F_Fmsy_true = numeric (n_years),HCR_status = character (n_years),stringsAsFactors = FALSE # История наблюдений для оценки <- numeric ()<- numeric ()<- NULL <- NULL # Основной цикл симуляции по годам for (t in 1 : n_years) {# 1. Генерация наблюдений (индекс биомассы) <- generate_index_observation ($ q, $ observation_error_cv<- c (index_history, index_obs)# 2. Оценка запаса (периодически или в первый год) if (t == 1 || (t - 1 ) %% settings$ assessment_interval == 0 ) {<- assess_stock_status ($ q,$ Bmsy,$ Fmsy,$ assessment_bias_cv# 3. Применение правила управления (HCR) <- HCR_function (last_assessment, previous_TAC)<- advice$ F_advice<- advice$ TAC<- TAC<- advice$ statuselse {# Используем прошлогодние рекомендации <- results$ F_advice[t-1 ]<- results$ TAC[t-1 ]<- results$ HCR_status[t-1 ]# 4. Реализация промысла с ошибкой implementation <- rlnorm (1 , meanlog = - settings$ implementation_error_cv^ 2 / 2 ,sdlog = settings$ implementation_error_cv<- F_advice * implementation_error<- min (F_realized, 2.0 ) # Ограничение максимальной F # 5. Расчет динамики популяции <- simulate_population_dynamics ($ process_error_cv# 6. Сохранение результатов текущего года $ B_true[t] <- B_sim$ B_estimated[t] <- ifelse (! is.null (last_assessment), $ B, NA )$ F_advice[t] <- F_advice$ F_realized[t] <- F_realized$ catch[t] <- pop_update$ catch_realized$ TAC[t] <- TAC$ index_obs[t] <- index_obs$ B_Bmsy_true[t] <- B_sim / params$ Bmsy$ F_Fmsy_true[t] <- F_realized / params$ Fmsy$ HCR_status[t] <- HCR_status# 7. Обновление состояния для следующего года <- pop_update$ B_next<- c (catch_history, pop_update$ catch_realized)return (results)# ------------------- 7. ЗАПУСК MSE ДЛЯ ВСЕХ СТРАТЕГИЙ -------------------- cat (" \n ========== ЗАПУСК СИМУЛЯЦИЙ MSE ========== \n " )
========== ЗАПУСК СИМУЛЯЦИЙ MSE ==========
## 7.1 Подготовка параметров для симуляций <- list (r = r_true,K = K_true,Bmsy = Bmsy_true,Fmsy = Fmsy_true,B_initial = B_current,q = q_cpue<- list (n_years = n_years,assessment_interval = assessment_interval,process_error_cv = process_error_cv,observation_error_cv = observation_error_cv,assessment_bias_cv = assessment_bias_cv,implementation_error_cv = implementation_error_cv## 7.2 Список стратегий для сравнения <- list ("Fish at Fmsy" = HCR_Fmsy,"MSY Hockey-stick" = HCR_hockey_stick,"ICES Advice Rule" = HCR_ICES## 7.3 Запуск симуляций для каждой стратегии <- list ()for (strategy_name in names (strategies)) {cat (sprintf (" \n Запуск стратегии: %s \n " , strategy_name))cat ("Прогресс: " )# Запуск n_sim симуляций для текущей стратегии <- list ()# Прогресс-индикатор <- txtProgressBar (min = 0 , max = n_sim, style = 3 , width = 50 )for (i in 1 : n_sim) {<- run_mse_simulation (sim_id = i,HCR_function = strategies[[strategy_name]],params = sim_params,settings = sim_settings$ strategy <- strategy_name<- sim_resultsetTxtProgressBar (pb, i)close (pb)# Объединение результатов симуляций <- bind_rows (strategy_results)cat (sprintf (" \n ✓ Завершено %d симуляций для %s \n " , n_sim, strategy_name))
Запуск стратегии: Fish at Fmsy
Прогресс:
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|= | 3%
|
|== | 3%
|
|== | 4%
|
|== | 5%
|
|=== | 5%
|
|=== | 6%
|
|=== | 7%
|
|==== | 7%
|
|==== | 8%
|
|==== | 9%
|
|===== | 9%
|
|===== | 10%
|
|===== | 11%
|
|====== | 11%
|
|====== | 12%
|
|====== | 13%
|
|======= | 13%
|
|======= | 14%
|
|======= | 15%
|
|======== | 15%
|
|======== | 16%
|
|======== | 17%
|
|========= | 17%
|
|========= | 18%
|
|========= | 19%
|
|========== | 19%
|
|========== | 20%
|
|========== | 21%
|
|=========== | 21%
|
|=========== | 22%
|
|=========== | 23%
|
|============ | 23%
|
|============ | 24%
|
|============ | 25%
|
|============= | 25%
|
|============= | 26%
|
|============= | 27%
|
|============== | 27%
|
|============== | 28%
|
|============== | 29%
|
|=============== | 29%
|
|=============== | 30%
|
|=============== | 31%
|
|================ | 31%
|
|================ | 32%
|
|================ | 33%
|
|================= | 33%
|
|================= | 34%
|
|================= | 35%
|
|================== | 35%
|
|================== | 36%
|
|================== | 37%
|
|=================== | 37%
|
|=================== | 38%
|
|=================== | 39%
|
|==================== | 39%
|
|==================== | 40%
|
|==================== | 41%
|
|===================== | 41%
|
|===================== | 42%
|
|===================== | 43%
|
|====================== | 43%
|
|====================== | 44%
|
|====================== | 45%
|
|======================= | 45%
|
|======================= | 46%
|
|======================= | 47%
|
|======================== | 47%
|
|======================== | 48%
|
|======================== | 49%
|
|========================= | 49%
|
|========================= | 50%
|
|========================= | 51%
|
|========================== | 51%
|
|========================== | 52%
|
|========================== | 53%
|
|=========================== | 53%
|
|=========================== | 54%
|
|=========================== | 55%
|
|============================ | 55%
|
|============================ | 56%
|
|============================ | 57%
|
|============================= | 57%
|
|============================= | 58%
|
|============================= | 59%
|
|============================== | 59%
|
|============================== | 60%
|
|============================== | 61%
|
|=============================== | 61%
|
|=============================== | 62%
|
|=============================== | 63%
|
|================================ | 63%
|
|================================ | 64%
|
|================================ | 65%
|
|================================= | 65%
|
|================================= | 66%
|
|================================= | 67%
|
|================================== | 67%
|
|================================== | 68%
|
|================================== | 69%
|
|=================================== | 69%
|
|=================================== | 70%
|
|=================================== | 71%
|
|==================================== | 71%
|
|==================================== | 72%
|
|==================================== | 73%
|
|===================================== | 73%
|
|===================================== | 74%
|
|===================================== | 75%
|
|====================================== | 75%
|
|====================================== | 76%
|
|====================================== | 77%
|
|======================================= | 77%
|
|======================================= | 78%
|
|======================================= | 79%
|
|======================================== | 79%
|
|======================================== | 80%
|
|======================================== | 81%
|
|========================================= | 81%
|
|========================================= | 82%
|
|========================================= | 83%
|
|========================================== | 83%
|
|========================================== | 84%
|
|========================================== | 85%
|
|=========================================== | 85%
|
|=========================================== | 86%
|
|=========================================== | 87%
|
|============================================ | 87%
|
|============================================ | 88%
|
|============================================ | 89%
|
|============================================= | 89%
|
|============================================= | 90%
|
|============================================= | 91%
|
|============================================== | 91%
|
|============================================== | 92%
|
|============================================== | 93%
|
|=============================================== | 93%
|
|=============================================== | 94%
|
|=============================================== | 95%
|
|================================================ | 95%
|
|================================================ | 96%
|
|================================================ | 97%
|
|================================================= | 97%
|
|================================================= | 98%
|
|================================================= | 99%
|
|==================================================| 99%
|
|==================================================| 100%
<U+2713> Завершено 500 симуляций для Fish at Fmsy
Запуск стратегии: MSY Hockey-stick
Прогресс:
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|= | 3%
|
|== | 3%
|
|== | 4%
|
|== | 5%
|
|=== | 5%
|
|=== | 6%
|
|=== | 7%
|
|==== | 7%
|
|==== | 8%
|
|==== | 9%
|
|===== | 9%
|
|===== | 10%
|
|===== | 11%
|
|====== | 11%
|
|====== | 12%
|
|====== | 13%
|
|======= | 13%
|
|======= | 14%
|
|======= | 15%
|
|======== | 15%
|
|======== | 16%
|
|======== | 17%
|
|========= | 17%
|
|========= | 18%
|
|========= | 19%
|
|========== | 19%
|
|========== | 20%
|
|========== | 21%
|
|=========== | 21%
|
|=========== | 22%
|
|=========== | 23%
|
|============ | 23%
|
|============ | 24%
|
|============ | 25%
|
|============= | 25%
|
|============= | 26%
|
|============= | 27%
|
|============== | 27%
|
|============== | 28%
|
|============== | 29%
|
|=============== | 29%
|
|=============== | 30%
|
|=============== | 31%
|
|================ | 31%
|
|================ | 32%
|
|================ | 33%
|
|================= | 33%
|
|================= | 34%
|
|================= | 35%
|
|================== | 35%
|
|================== | 36%
|
|================== | 37%
|
|=================== | 37%
|
|=================== | 38%
|
|=================== | 39%
|
|==================== | 39%
|
|==================== | 40%
|
|==================== | 41%
|
|===================== | 41%
|
|===================== | 42%
|
|===================== | 43%
|
|====================== | 43%
|
|====================== | 44%
|
|====================== | 45%
|
|======================= | 45%
|
|======================= | 46%
|
|======================= | 47%
|
|======================== | 47%
|
|======================== | 48%
|
|======================== | 49%
|
|========================= | 49%
|
|========================= | 50%
|
|========================= | 51%
|
|========================== | 51%
|
|========================== | 52%
|
|========================== | 53%
|
|=========================== | 53%
|
|=========================== | 54%
|
|=========================== | 55%
|
|============================ | 55%
|
|============================ | 56%
|
|============================ | 57%
|
|============================= | 57%
|
|============================= | 58%
|
|============================= | 59%
|
|============================== | 59%
|
|============================== | 60%
|
|============================== | 61%
|
|=============================== | 61%
|
|=============================== | 62%
|
|=============================== | 63%
|
|================================ | 63%
|
|================================ | 64%
|
|================================ | 65%
|
|================================= | 65%
|
|================================= | 66%
|
|================================= | 67%
|
|================================== | 67%
|
|================================== | 68%
|
|================================== | 69%
|
|=================================== | 69%
|
|=================================== | 70%
|
|=================================== | 71%
|
|==================================== | 71%
|
|==================================== | 72%
|
|==================================== | 73%
|
|===================================== | 73%
|
|===================================== | 74%
|
|===================================== | 75%
|
|====================================== | 75%
|
|====================================== | 76%
|
|====================================== | 77%
|
|======================================= | 77%
|
|======================================= | 78%
|
|======================================= | 79%
|
|======================================== | 79%
|
|======================================== | 80%
|
|======================================== | 81%
|
|========================================= | 81%
|
|========================================= | 82%
|
|========================================= | 83%
|
|========================================== | 83%
|
|========================================== | 84%
|
|========================================== | 85%
|
|=========================================== | 85%
|
|=========================================== | 86%
|
|=========================================== | 87%
|
|============================================ | 87%
|
|============================================ | 88%
|
|============================================ | 89%
|
|============================================= | 89%
|
|============================================= | 90%
|
|============================================= | 91%
|
|============================================== | 91%
|
|============================================== | 92%
|
|============================================== | 93%
|
|=============================================== | 93%
|
|=============================================== | 94%
|
|=============================================== | 95%
|
|================================================ | 95%
|
|================================================ | 96%
|
|================================================ | 97%
|
|================================================= | 97%
|
|================================================= | 98%
|
|================================================= | 99%
|
|==================================================| 99%
|
|==================================================| 100%
<U+2713> Завершено 500 симуляций для MSY Hockey-stick
Запуск стратегии: ICES Advice Rule
Прогресс:
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|= | 3%
|
|== | 3%
|
|== | 4%
|
|== | 5%
|
|=== | 5%
|
|=== | 6%
|
|=== | 7%
|
|==== | 7%
|
|==== | 8%
|
|==== | 9%
|
|===== | 9%
|
|===== | 10%
|
|===== | 11%
|
|====== | 11%
|
|====== | 12%
|
|====== | 13%
|
|======= | 13%
|
|======= | 14%
|
|======= | 15%
|
|======== | 15%
|
|======== | 16%
|
|======== | 17%
|
|========= | 17%
|
|========= | 18%
|
|========= | 19%
|
|========== | 19%
|
|========== | 20%
|
|========== | 21%
|
|=========== | 21%
|
|=========== | 22%
|
|=========== | 23%
|
|============ | 23%
|
|============ | 24%
|
|============ | 25%
|
|============= | 25%
|
|============= | 26%
|
|============= | 27%
|
|============== | 27%
|
|============== | 28%
|
|============== | 29%
|
|=============== | 29%
|
|=============== | 30%
|
|=============== | 31%
|
|================ | 31%
|
|================ | 32%
|
|================ | 33%
|
|================= | 33%
|
|================= | 34%
|
|================= | 35%
|
|================== | 35%
|
|================== | 36%
|
|================== | 37%
|
|=================== | 37%
|
|=================== | 38%
|
|=================== | 39%
|
|==================== | 39%
|
|==================== | 40%
|
|==================== | 41%
|
|===================== | 41%
|
|===================== | 42%
|
|===================== | 43%
|
|====================== | 43%
|
|====================== | 44%
|
|====================== | 45%
|
|======================= | 45%
|
|======================= | 46%
|
|======================= | 47%
|
|======================== | 47%
|
|======================== | 48%
|
|======================== | 49%
|
|========================= | 49%
|
|========================= | 50%
|
|========================= | 51%
|
|========================== | 51%
|
|========================== | 52%
|
|========================== | 53%
|
|=========================== | 53%
|
|=========================== | 54%
|
|=========================== | 55%
|
|============================ | 55%
|
|============================ | 56%
|
|============================ | 57%
|
|============================= | 57%
|
|============================= | 58%
|
|============================= | 59%
|
|============================== | 59%
|
|============================== | 60%
|
|============================== | 61%
|
|=============================== | 61%
|
|=============================== | 62%
|
|=============================== | 63%
|
|================================ | 63%
|
|================================ | 64%
|
|================================ | 65%
|
|================================= | 65%
|
|================================= | 66%
|
|================================= | 67%
|
|================================== | 67%
|
|================================== | 68%
|
|================================== | 69%
|
|=================================== | 69%
|
|=================================== | 70%
|
|=================================== | 71%
|
|==================================== | 71%
|
|==================================== | 72%
|
|==================================== | 73%
|
|===================================== | 73%
|
|===================================== | 74%
|
|===================================== | 75%
|
|====================================== | 75%
|
|====================================== | 76%
|
|====================================== | 77%
|
|======================================= | 77%
|
|======================================= | 78%
|
|======================================= | 79%
|
|======================================== | 79%
|
|======================================== | 80%
|
|======================================== | 81%
|
|========================================= | 81%
|
|========================================= | 82%
|
|========================================= | 83%
|
|========================================== | 83%
|
|========================================== | 84%
|
|========================================== | 85%
|
|=========================================== | 85%
|
|=========================================== | 86%
|
|=========================================== | 87%
|
|============================================ | 87%
|
|============================================ | 88%
|
|============================================ | 89%
|
|============================================= | 89%
|
|============================================= | 90%
|
|============================================= | 91%
|
|============================================== | 91%
|
|============================================== | 92%
|
|============================================== | 93%
|
|=============================================== | 93%
|
|=============================================== | 94%
|
|=============================================== | 95%
|
|================================================ | 95%
|
|================================================ | 96%
|
|================================================ | 97%
|
|================================================= | 97%
|
|================================================= | 98%
|
|================================================= | 99%
|
|==================================================| 99%
|
|==================================================| 100%
<U+2713> Завершено 500 симуляций для ICES Advice Rule
## 7.4 Объединение всех результатов в один датафрейм <- bind_rows (all_results)cat (" \n ✓ Всего выполнено:" , n_sim * length (strategies), "симуляций \n " )
<U+2713> Всего выполнено: 1500 симуляций
# ------------------- 8. РАСЧЕТ РИСК-МЕТРИК -------------------- cat (" \n ========== РАСЧЕТ РИСК-МЕТРИК ========== \n " )
========== РАСЧЕТ РИСК-МЕТРИК ==========
## 8.1 Функция для расчета комплексных метрик производительности <- function (results) {<- results %>% group_by (strategy) %>% summarise (# === БИОЛОГИЧЕСКИЕ МЕТРИКИ === # Вероятность перелова (F > Fmsy) prob_overfishing = mean (F_Fmsy_true > 1 , na.rm = TRUE ),# Вероятность истощения запаса (B < 0.5*Bmsy) prob_overfished = mean (B_Bmsy_true < 0.5 , na.rm = TRUE ),# Вероятность коллапса (B < 0.2*Bmsy) prob_collapsed = mean (B_Bmsy_true < 0.2 , na.rm = TRUE ),# Вероятность нахождения в "зеленой зоне" Кобе prob_green_zone = mean (B_Bmsy_true > 1 & F_Fmsy_true < 1 , na.rm = TRUE ),# Средние показатели за последние 5 лет mean_B_Bmsy_final = mean (B_Bmsy_true[year > (max (year) - 5 )], na.rm = TRUE ),mean_F_Fmsy_final = mean (F_Fmsy_true[year > (max (year) - 5 )], na.rm = TRUE ),# === ЭКОНОМИЧЕСКИЕ МЕТРИКИ === # Средний вылов за весь период mean_catch = mean (catch, na.rm = TRUE ),# Суммарный вылов за период total_catch = sum (catch) / n_distinct (sim_id),# Стабильность вылова (обратная величина коэффициента вариации) catch_stability = {<- aggregate (catch, by = list (year), mean)$ x1 - (sd (annual_catch, na.rm = TRUE ) / mean (annual_catch, na.rm = TRUE ))# Средняя межгодовая изменчивость вылова (AAV) catch_aav = {<- aggregate (catch, by = list (year), mean)$ xmean (abs (diff (annual_catch)) / annual_catch[- length (annual_catch)], na.rm = TRUE )# === УПРАВЛЕНЧЕСКИЕ МЕТРИКИ === # Вероятность закрытия промысла prob_closure = mean (TAC == 0 , na.rm = TRUE ),# Частота изменения рекомендаций advice_changes = {<- aggregate (F_advice, by = list (sim_id), function (x) sum (diff (x) != 0 ))$ xmean (changes) / (max (year) - 1 ).groups = "drop" return (metrics)## 8.2 Расчет метрик для каждой стратегии <- calculate_performance_metrics (mse_results)## 8.3 Добавление комплексных оценок <- performance_metrics %>% mutate (# Биологический индекс (0-1, где 1 - лучше) bio_score = 0.4 * (1 - prob_overfished) + 0.3 * (1 - prob_overfishing) + 0.3 * prob_green_zone,# Экономический индекс (0-1, где 1 - лучше) econ_score = 0.5 * (mean_catch / max (mean_catch)) + 0.3 * catch_stability + 0.2 * (1 - catch_aav),# Общий индекс производительности overall_score = 0.6 * bio_score + 0.4 * econ_score%>% arrange (desc (overall_score))## 8.4 Вывод таблицы метрик cat (" \n --- ТАБЛИЦА РИСК-МЕТРИК --- \n\n " )
--- ТАБЛИЦА РИСК-МЕТРИК ---
# Форматированный вывод ключевых метрик for (strat in unique (performance_metrics$ strategy)) {<- performance_metrics[performance_metrics$ strategy == strat, ]cat (sprintf ("СТРАТЕГИЯ: %s \n " , strat))cat (strrep ("-" , 40 ), " \n " )cat (sprintf ("Биологические риски: \n " ))cat (sprintf (" P(перелов): %.1f%% \n " , metrics_row$ prob_overfishing * 100 ))cat (sprintf (" P(истощение): %.1f%% \n " , metrics_row$ prob_overfished * 100 ))cat (sprintf (" P(коллапс): %.1f%% \n " , metrics_row$ prob_collapsed * 100 ))cat (sprintf (" P(зеленая зона): %.1f%% \n " , metrics_row$ prob_green_zone * 100 ))cat (sprintf ("Экономические показатели: \n " ))cat (sprintf (" Средний вылов: %.1f тыс. т \n " , metrics_row$ mean_catch))cat (sprintf (" Стабильность: %.2f \n " , metrics_row$ catch_stability))cat (sprintf (" Межгодовая изменчивость: %.1f%% \n " , metrics_row$ catch_aav * 100 ))cat (sprintf ("Итоговые оценки: \n " ))cat (sprintf (" Биологический индекс: %.2f \n " , metrics_row$ bio_score))cat (sprintf (" Экономический индекс: %.2f \n " , metrics_row$ econ_score))cat (sprintf (" ОБЩИЙ ИНДЕКС: %.2f \n " , metrics_row$ overall_score))cat (" \n " )
СТРАТЕГИЯ: ICES Advice Rule
----------------------------------------
Биологические риски:
P(перелов): 8.3%
P(истощение): 0.2%
P(коллапс): 0.0%
P(зеленая зона): 57.6%
Экономические показатели:
Средний вылов: 10.3 тыс. т
Стабильность: 0.95
Межгодовая изменчивость: 1.3%
Итоговые оценки:
Биологический индекс: 0.85
Экономический индекс: 0.97
ОБЩИЙ ИНДЕКС: 0.90
СТРАТЕГИЯ: MSY Hockey-stick
----------------------------------------
Биологические риски:
P(перелов): 46.3%
P(истощение): 3.7%
P(коллапс): 0.0%
P(зеленая зона): 20.2%
Экономические показатели:
Средний вылов: 10.5 тыс. т
Стабильность: 0.92
Межгодовая изменчивость: 1.1%
Итоговые оценки:
Биологический индекс: 0.61
Экономический индекс: 0.97
ОБЩИЙ ИНДЕКС: 0.75
СТРАТЕГИЯ: Fish at Fmsy
----------------------------------------
Биологические риски:
P(перелов): 48.5%
P(истощение): 4.7%
P(коллапс): 0.0%
P(зеленая зона): 20.0%
Экономические показатели:
Средний вылов: 10.5 тыс. т
Стабильность: 0.92
Межгодовая изменчивость: 1.0%
Итоговые оценки:
Биологический индекс: 0.60
Экономический индекс: 0.97
ОБЩИЙ ИНДЕКС: 0.75
# ------------------- 9. ВИЗУАЛИЗАЦИЯ РЕЗУЛЬТАТОВ -------------------- cat ("========== СОЗДАНИЕ ГРАФИКОВ ========== \n " )
========== СОЗДАНИЕ ГРАФИКОВ ==========
## 9.1 Настройка темы для графиков <- theme_minimal () + theme (plot.title = element_text (size = 14 , face = "bold" ),plot.subtitle = element_text (size = 11 ),strip.text = element_text (size = 11 , face = "bold" ),legend.position = "bottom" ,panel.grid.minor = element_blank ()## 9.2 Подготовка сводных данных для визуализации <- mse_results %>% group_by (strategy, year) %>% summarise (# Квантили для биомассы B_median = median (B_true, na.rm = TRUE ),B_q25 = quantile (B_true, 0.25 , na.rm = TRUE ),B_q75 = quantile (B_true, 0.75 , na.rm = TRUE ),B_q05 = quantile (B_true, 0.05 , na.rm = TRUE ),B_q95 = quantile (B_true, 0.95 , na.rm = TRUE ),# Квантили для B/Bmsy B_Bmsy_median = median (B_Bmsy_true, na.rm = TRUE ),B_Bmsy_q25 = quantile (B_Bmsy_true, 0.25 , na.rm = TRUE ),B_Bmsy_q75 = quantile (B_Bmsy_true, 0.75 , na.rm = TRUE ),# Квантили для F/Fmsy F_Fmsy_median = median (F_Fmsy_true, na.rm = TRUE ),F_Fmsy_q25 = quantile (F_Fmsy_true, 0.25 , na.rm = TRUE ),F_Fmsy_q75 = quantile (F_Fmsy_true, 0.75 , na.rm = TRUE ),# Квантили для вылова catch_median = median (catch, na.rm = TRUE ),catch_q25 = quantile (catch, 0.25 , na.rm = TRUE ),catch_q75 = quantile (catch, 0.75 , na.rm = TRUE ),.groups = "drop" ## 9.3 График 1: Траектории B/Bmsy <- ggplot (summary_data, aes (x = year, color = strategy, fill = strategy)) + geom_ribbon (aes (ymin = B_Bmsy_q25, ymax = B_Bmsy_q75), alpha = 0.3 , color = NA ) + geom_line (aes (y = B_Bmsy_median), size = 1.5 ) + geom_hline (yintercept = 1 , linetype = "dashed" , color = "black" , alpha = 0.5 ) + geom_hline (yintercept = 0.5 , linetype = "dotted" , color = "red" , alpha = 0.5 ) + scale_color_manual (values = c ("Fish at Fmsy" = "#E41A1C" , "MSY Hockey-stick" = "#377EB8" ,"ICES Advice Rule" = "#4DAF4A" )) + scale_fill_manual (values = c ("Fish at Fmsy" = "#E41A1C" , "MSY Hockey-stick" = "#377EB8" ,"ICES Advice Rule" = "#4DAF4A" )) + labs (title = "Динамика относительной биомассы (B/Bmsy)" ,subtitle = "Медиана и межквартильный размах" ,x = "Год прогноза" , y = "B/Bmsy" ,color = "Стратегия" ,fill = "Стратегия" ) + + annotate ("text" , x = max (summary_data$ year), y = 1.02 , label = "Bmsy" , hjust = 1 , vjust = 0 , size = 3 ) + annotate ("text" , x = max (summary_data$ year), y = 0.52 , label = "0.5 Bmsy" , hjust = 1 , vjust = 0 , size = 3 , color = "red" )
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
i Please use `linewidth` instead.
## 9.4 График 2: Траектории F/Fmsy <- ggplot (summary_data, aes (x = year, color = strategy, fill = strategy)) + geom_ribbon (aes (ymin = F_Fmsy_q25, ymax = F_Fmsy_q75), alpha = 0.3 , color = NA ) + geom_line (aes (y = F_Fmsy_median), size = 1.5 ) + geom_hline (yintercept = 1 , linetype = "dashed" , color = "black" , alpha = 0.5 ) + scale_color_manual (values = c ("Fish at Fmsy" = "#E41A1C" , "MSY Hockey-stick" = "#377EB8" ,"ICES Advice Rule" = "#4DAF4A" )) + scale_fill_manual (values = c ("Fish at Fmsy" = "#E41A1C" , "MSY Hockey-stick" = "#377EB8" ,"ICES Advice Rule" = "#4DAF4A" )) + labs (title = "Динамика промысловой смертности (F/Fmsy)" ,subtitle = "Медиана и межквартильный размах" ,x = "Год прогноза" , y = "F/Fmsy" ,color = "Стратегия" ,fill = "Стратегия" ) + + coord_cartesian (ylim = c (0 , 2 )) + annotate ("text" , x = max (summary_data$ year), y = 1.02 , label = "Fmsy" , hjust = 1 , vjust = 0 , size = 3 )## 9.5 График 3: Динамика вылова <- ggplot (summary_data, aes (x = year, color = strategy, fill = strategy)) + geom_ribbon (aes (ymin = catch_q25, ymax = catch_q75), alpha = 0.3 , color = NA ) + geom_line (aes (y = catch_median), size = 1.5 ) + scale_color_manual (values = c ("Fish at Fmsy" = "#E41A1C" , "MSY Hockey-stick" = "#377EB8" ,"ICES Advice Rule" = "#4DAF4A" )) + scale_fill_manual (values = c ("Fish at Fmsy" = "#E41A1C" , "MSY Hockey-stick" = "#377EB8" ,"ICES Advice Rule" = "#4DAF4A" )) + labs (title = "Динамика вылова" ,subtitle = "Медиана и межквартильный размах" ,x = "Год прогноза" , y = "Вылов (тыс. т)" ,color = "Стратегия" ,fill = "Стратегия" ) + ## 9.6 График 4: Фазовая диаграмма Кобе для последнего года <- mse_results %>% filter (year == max (year)) %>% group_by (strategy) %>% sample_n (min (100 , n ())) # Максимум 100 точек для читаемости <- ggplot (kobe_data, aes (x = B_Bmsy_true, y = F_Fmsy_true)) + # Зоны Кобе annotate ("rect" , xmin = 0 , xmax = 1 , ymin = 1 , ymax = 3 ,fill = "red" , alpha = 0.15 ) + # Красная зона annotate ("rect" , xmin = 1 , xmax = 3 , ymin = 1 , ymax = 3 ,fill = "#FFA500" , alpha = 0.15 ) + # Оранжевая зона annotate ("rect" , xmin = 0 , xmax = 1 , ymin = 0 , ymax = 1 ,fill = "#FFFF00" , alpha = 0.15 ) + # Желтая зона annotate ("rect" , xmin = 1 , xmax = 3 , ymin = 0 , ymax = 1 ,fill = "green" , alpha = 0.15 ) + # Зеленая зона # Точки симуляций geom_point (aes (color = strategy), alpha = 0.5 , size = 1.5 ) + # Референсные линии geom_vline (xintercept = 1 , linetype = "solid" , size = 0.8 ) + geom_hline (yintercept = 1 , linetype = "solid" , size = 0.8 ) + # Настройки facet_wrap (~ strategy, ncol = 3 ) + scale_color_manual (values = c ("Fish at Fmsy" = "#E41A1C" , "MSY Hockey-stick" = "#377EB8" ,"ICES Advice Rule" = "#4DAF4A" )) + labs (title = "Фазовая диаграмма Кобе (год 20)" ,subtitle = "Распределение конечных состояний" ,x = "B/Bmsy" , y = "F/Fmsy" ) + + theme (legend.position = "none" ) + coord_cartesian (xlim = c (0 , 2 ), ylim = c (0 , 2 ))## 9.7 График 5: Сравнение риск-метрик (барплот) <- performance_metrics %>% select (strategy, prob_overfishing, prob_overfished, %>% pivot_longer (cols = - strategy, names_to = "metric" , values_to = "probability" )# Переименование метрик для графика <- c (prob_overfishing = "Перелов \n (F > Fmsy)" ,prob_overfished = "Истощение \n (B < 0.5 Bmsy)" ,prob_collapsed = "Коллапс \n (B < 0.2 Bmsy)" ,prob_green_zone = "Зеленая зона \n (B > Bmsy, F < Fmsy)" <- ggplot (metrics_long, aes (x = metric, y = probability, fill = strategy)) + geom_bar (stat = "identity" , position = position_dodge (width = 0.8 ), width = 0.7 ) + geom_text (aes (label = sprintf ("%.1f%%" , probability * 100 )),position = position_dodge (width = 0.8 ), vjust = - 0.5 , size = 3 ) + scale_fill_manual (values = c ("Fish at Fmsy" = "#E41A1C" , "MSY Hockey-stick" = "#377EB8" ,"ICES Advice Rule" = "#4DAF4A" )) + scale_x_discrete (labels = metric_labels) + scale_y_continuous (labels = scales:: percent, limits = c (0 , max (metrics_long$ probability) * 1.1 )) + labs (title = "Сравнение вероятностей риска" ,subtitle = "По всему периоду симуляции" ,x = "" , y = "Вероятность" ,fill = "Стратегия" ) + ## 9.8 График 6: Компромиссы (trade-offs) <- ggplot (performance_metrics, aes (x = mean_catch, y = 1 - prob_overfished)) + geom_point (aes (color = strategy, size = catch_stability), alpha = 0.7 ) + geom_text (aes (label = strategy), vjust = - 1.5 , size = 3.5 ) + scale_color_manual (values = c ("Fish at Fmsy" = "#E41A1C" , "MSY Hockey-stick" = "#377EB8" ,"ICES Advice Rule" = "#4DAF4A" )) + scale_size_continuous (range = c (8 , 15 ), name = "Стабильность \n вылова" ) + labs (title = "Анализ компромиссов между целями управления" ,subtitle = "Размер точки = стабильность вылова" ,x = "Средний вылов (тыс. т)" , y = "Вероятность избежать истощения" ,color = "Стратегия" ) + + theme (legend.position = "right" )# ------------------- 10. ДЕТАЛЬНЫЙ АНАЛИЗ -------------------- cat (" \n ========== ДЕТАЛЬНЫЙ АНАЛИЗ РЕЗУЛЬТАТОВ ========== \n " )
========== ДЕТАЛЬНЫЙ АНАЛИЗ РЕЗУЛЬТАТОВ ==========
## 10.1 Анализ времени восстановления cat (" \n --- Анализ восстановления запаса --- \n " )
--- Анализ восстановления запаса ---
if (B_current/ Bmsy_true < 1 ) {<- mse_results %>% group_by (strategy, sim_id) %>% summarise (recovery_year = which (B_Bmsy_true > 1 )[1 ],recovered = ! is.na (recovery_year),.groups = "drop" %>% group_by (strategy) %>% summarise (prob_recovery = mean (recovered, na.rm = TRUE ),median_recovery_time = median (recovery_year, na.rm = TRUE ),q25_recovery = quantile (recovery_year, 0.25 , na.rm = TRUE ),q75_recovery = quantile (recovery_year, 0.75 , na.rm = TRUE ),.groups = "drop" for (i in 1 : nrow (recovery_analysis)) {cat (sprintf ("%s: \n " , recovery_analysis$ strategy[i]))cat (sprintf (" Вероятность восстановления: %.1f%% \n " , $ prob_recovery[i] * 100 ))if (! is.na (recovery_analysis$ median_recovery_time[i])) {cat (sprintf (" Медианное время восстановления: %.0f лет [%.0f-%.0f] \n " ,$ median_recovery_time[i],$ q25_recovery[i],$ q75_recovery[i]))else {cat ("Запас находится на уровне B/Bmsy >= 1. Восстановление не требуется. \n " )cat (sprintf ("Текущее состояние: B/Bmsy = %.2f \n " , B_current / Bmsy_true))
Запас находится на уровне B/Bmsy >= 1. Восстановление не требуется.
Текущее состояние: B/Bmsy = 1.34
## 10.2 Анализ частоты закрытия промысла <- mse_results %>% group_by (strategy) %>% summarise (total_closures = sum (TAC == 0 ),closure_rate = mean (TAC == 0 ),avg_closure_duration = {<- rle (TAC == 0 )if (any (rle_results$ values)) {mean (rle_results$ lengths[rle_results$ values])else {0 .groups = "drop" cat (" \n --- Анализ закрытия промысла --- \n " )
--- Анализ закрытия промысла ---
for (i in 1 : nrow (closure_analysis)) {cat (sprintf ("%s: \n " , closure_analysis$ strategy[i]))cat (sprintf (" Частота закрытия: %.1f%% \n " , closure_analysis$ closure_rate[i] * 100 ))if (closure_analysis$ avg_closure_duration[i] > 0 ) {cat (sprintf (" Средняя продолжительность закрытия: %.1f лет \n " ,$ avg_closure_duration[i]))
Fish at Fmsy:
Частота закрытия: 0.0%
ICES Advice Rule:
Частота закрытия: 1.1%
Средняя продолжительность закрытия: 2.3 лет
MSY Hockey-stick:
Частота закрытия: 0.0%
Средняя продолжительность закрытия: 2.0 лет
# ------------------- 11. СОХРАНЕНИЕ РЕЗУЛЬТАТОВ -------------------- cat (" \n ========== СОХРАНЕНИЕ РЕЗУЛЬТАТОВ ========== \n " )
========== СОХРАНЕНИЕ РЕЗУЛЬТАТОВ ==========
## 11.2 Сохранение таблиц write.csv (performance_metrics, "MSE_performance_metrics.csv" , row.names = FALSE )write.csv (summary_data, "MSE_summary_by_year.csv" , row.names = FALSE )cat ("✓ Таблицы метрик сохранены \n " )
<U+2713> Таблицы метрик сохранены
## 11.3 Сохранение полных результатов saveRDS (mse_results, "MSE_full_results.rds" )cat ("✓ Полные результаты симуляций сохранены: 'MSE_full_results.rds' \n " )
<U+2713> Полные результаты симуляций сохранены: 'MSE_full_results.rds'
## 11.4 Создание текстового отчета sink ("MSE_report.txt" )cat ("ОТЧЕТ ПО MANAGEMENT STRATEGY EVALUATION \n " )cat (strrep ("=" , 60 ), " \n\n " )cat ("Дата анализа:" , format (Sys.Date (), "%d.%m.%Y" ), " \n " )cat ("Версия SPiCT:" , as.character (packageVersion ("spict" )), " \n\n " )cat ("ПАРАМЕТРЫ СИМУЛЯЦИИ: \n " )cat (sprintf (" Количество симуляций: %d \n " , n_sim))cat (sprintf (" Период прогноза: %d лет \n " , n_years))cat (sprintf (" Интервал оценки: каждые %d года \n " , assessment_interval))cat (sprintf (" CV процессной ошибки: %.0f%% \n " , process_error_cv * 100 ))cat (sprintf (" CV ошибки наблюдения: %.0f%% \n " , observation_error_cv * 100 ))cat (sprintf (" CV смещения оценок: %.0f%% \n " , assessment_bias_cv * 100 ))cat (sprintf (" CV ошибки реализации: %.0f%% \n\n " , implementation_error_cv * 100 ))cat ("ИСХОДНОЕ СОСТОЯНИЕ ЗАПАСА: \n " )cat (sprintf (" B/Bmsy: %.2f \n " , B_current/ Bmsy_true))cat (sprintf (" F/Fmsy: %.2f \n\n " , F_current/ Fmsy_true))cat ("РЕЗУЛЬТАТЫ ПО СТРАТЕГИЯМ: \n\n " )for (i in 1 : nrow (performance_metrics)) {cat (sprintf ("СТРАТЕГИЯ: %s \n " , performance_metrics$ strategy[i]))cat (strrep ("-" , 40 ), " \n " )cat (sprintf (" Риск перелова: %.1f%% \n " , performance_metrics$ prob_overfishing[i] * 100 ))cat (sprintf (" Риск истощения: %.1f%% \n " , performance_metrics$ prob_overfished[i] * 100 ))cat (sprintf (" Вероятность устойчивого состояния: %.1f%% \n " , $ prob_green_zone[i] * 100 ))cat (sprintf (" Средний вылов: %.1f тыс. т \n " , performance_metrics$ mean_catch[i]))cat (sprintf (" Стабильность вылова: %.2f \n " , performance_metrics$ catch_stability[i]))cat (sprintf (" Общий индекс производительности: %.3f \n\n " , $ overall_score[i]))cat ("РЕКОМЕНДАЦИИ: \n " )<- performance_metrics$ strategy[1 ]cat (sprintf ("Оптимальная стратегия: %s \n " , best_strategy))cat (" \n Обоснование: \n " )if (best_strategy == "ICES Advice Rule" ) {cat ("- Предосторожный подход минимизирует биологические риски \n " )cat ("- Ограничение межгодовых изменений обеспечивает стабильность для промышленности \n " )cat ("- Баланс между сохранением запаса и экономической эффективностью \n " )else if (best_strategy == "MSY Hockey-stick" ) {cat ("- Адаптивное управление в зависимости от состояния запаса \n " )cat ("- Защита запаса при низкой биомассе \n " )cat ("- Возможность полного использования при хорошем состоянии \n " )else {cat ("- Простота применения \n " )cat ("- Максимизация долгосрочного вылова \n " )cat ("- Требует точной оценки Fmsy \n " )sink ()cat ("✓ Текстовый отчет сохранен: 'MSE_report.txt' \n " )
<U+2713> Текстовый отчет сохранен: 'MSE_report.txt'
cat (" \n " , strrep ("=" , 60 ), " \n " )
============================================================
cat (" MSE АНАЛИЗ УСПЕШНО ЗАВЕРШЕН \n " )
MSE АНАЛИЗ УСПЕШНО ЗАВЕРШЕН
cat (strrep ("=" , 60 ), " \n " )
============================================================
Результаты Management Strategy Evaluation: Сравнение стратегий управления Результаты выполнения третьего скрипта представляют собой комплексный анализ эффективности трех ключевых стратегий управления в условиях неопределенности. После 500 симуляций для каждой стратегии на 100-летнем горизонте получены следующие результаты.
Стратегия “Fish at Fmsy” показала наихудшие биологические показатели. Вероятность перелова составила 48.5%, вероятность истощения запаса (B < 0.5Bmsy) — 4.7%. Несмотря на наибольший средний вылов (10.5 тыс. тонн), эта стратегия демонстрирует неприемлемо высокие риски. Только 20% времени запас находился в “зеленой зоне” — состоянии, когда биомасса превышает Bmsy при промысловой смертности ниже Fmsy.
Стратегия “MSY Hockey-stick” показала умеренное улучшение. Вероятность перелова снизилась до 46.3%, вероятность истощения — до 3.7%. Однако показатель нахождения в зеленой зоне остался на низком уровне (20.2%). Средний вылов практически не изменился (10.5 тыс. тонн), что указывает на ограниченную эффективность данного подхода в конкретных условиях моделирования.
Стратегия “ICES Advice Rule” продемонстрировала наилучшие результаты. Вероятность перелова составила лишь 8.3%, вероятность истощения — 0.2%. Запас проводил 57.6% времени в зеленой зоне, что свидетельствует о высокой устойчивости. При этом средний вылов (10.3 тыс. тонн) практически не уступил другим стратегиям.
Анализ динамики показал, что адаптивные стратегии обеспечивают более стабильное состояние запаса в долгосрочной перспективе. К 20-му году прогноза стратегия ICES стабилизировала запас на уровне 1.0-1.2 B/Bmsy, в то время как другие стратегии демонстрировали значительные колебания.
Экономические показатели также свидетельствуют в пользу предосторожного подхода. Стабильность вылова при стратегии ICES составила 0.95 против 0.92 у других стратегий, при этом межгодовая изменчивость вылова была минимальной (1.3%).
Общий индекс производительности, учитывающий как биологические, так и экономические аспекты, составил 0.90 для ICES против 0.75 для обеих других стратегий. Это указывает на явное преимущество предосторожного подхода в условиях неопределенности.
Как мог бы заключить Довлатов, “идеальное управление — это когда и рыба цела, и рыбаки довольны, но в реальности приходится выбирать между риском и уверенностью”. Результаты MSE однозначно показывают, что стратегия ICES Advice Rule обеспечивает оптимальный баланс между сохранением ресурса и экономической эффективностью.
Важно отметить, что эти результаты справедливы для конкретных условий моделирования — начального состояния запаса (B/Bmsy = 1.34) и принятых уровней неопределенности. Изменение исходных параметров может повлиять на сравнительную эффективность стратегий, однако общее преимущество предосторожного подхода сохраняется.
