第五章：数学模型与控制理论

// SPDX-License-Identifier: MIT
pragma solidity ^0.8.0;

contract PIDController {
    // PID参数
    int256 public Kp = 1e16;      // 比例系数 (0.01)
    int256 public Ki = 1e15;      // 积分系数 (0.001)
    int256 public Kd = 5e15;      // 微分系数 (0.005)
    
    // 状态变量
    int256 public integral;        // 积分项
    int256 public lastError;       // 上次误差
    uint256 public lastTime;       // 上次更新时间
    
    // 目标价格和界限
    uint256 public constant TARGET_PRICE = 1e18;  // 1 USD
    int256 public constant MAX_OUTPUT = 5e16;     // 最大调整5%
    int256 public constant MIN_OUTPUT = -5e16;    // 最小调整-5%
    
    // 死区设置（减少频繁调整）
    uint256 public constant DEADBAND = 1e15;      // 0.1%死区
    
    // 计算PID输出
    function computePID(uint256 currentPrice) public returns (int256) {
        require(currentPrice > 0, "Invalid price");
        
        uint256 currentTime = block.timestamp;
        uint256 deltaTime = currentTime - lastTime;
        
        // 计算误差（可以为负）
        int256 error = int256(TARGET_PRICE) - int256(currentPrice);
        
        // 死区处理
        if (abs(error) < int256(DEADBAND)) {
            return 0;
        }
        
        // 计算PID各项
        int256 proportional = Kp * error / 1e18;
        
        // 积分项（带抗积分饱和）
        integral += error * int256(deltaTime);
        // 限制积分项防止积分饱和
        if (integral > MAX_OUTPUT * 1e18 / Ki) {
            integral = MAX_OUTPUT * 1e18 / Ki;
        } else if (integral < MIN_OUTPUT * 1e18 / Ki) {
            integral = MIN_OUTPUT * 1e18 / Ki;
        }
        int256 integralTerm = Ki * integral / 1e18;
        
        // 微分项（带滤波）
        int256 derivative = 0;
        if (deltaTime > 0) {
            derivative = Kd * (error - lastError) * 1e18 / int256(deltaTime);
        }
        
        // 计算总输出
        int256 output = proportional + integralTerm + derivative;
        
        // 输出限幅
        if (output > MAX_OUTPUT) {
            output = MAX_OUTPUT;
        } else if (output < MIN_OUTPUT) {
            output = MIN_OUTPUT;
        }
        
        // 更新状态
        lastError = error;
        lastTime = currentTime;
        
        emit PIDComputed(currentPrice, error, output);
        
        return output;
    }
    
    // 参数自适应调整
    function adaptiveParameterTuning(
        uint256 volatility,
        uint256 volume
    ) external {
        // 高波动性时增加Kd，减少Ki
        if (volatility > 2e16) { // 2%波动率
            Kd = Kd * 12 / 10;    // 增加20%
            Ki = Ki * 8 / 10;     // 减少20%
        }
        
        // 低交易量时减少所有增益
        if (volume < 1e6 * 1e18) { // 少于100万美元
            Kp = Kp * 9 / 10;
            Ki = Ki * 9 / 10;
            Kd = Kd * 9 / 10;
        }
    }
    
    // 辅助函数：计算绝对值
    function abs(int256 x) private pure returns (int256) {
        return x >= 0 ? x : -x;
    }
    
    event PIDComputed(uint256 price, int256 error, int256 output);
}

contract StabilityFeeController {
    PIDController public pidController;
    IPriceOracle public priceOracle;
    
    // 费率参数
    uint256 public baseRate = 2e16;        // 基础费率2%
    uint256 public currentRate;            // 当前费率
    uint256 public minRate = 0;            // 最低费率0%
    uint256 public maxRate = 20e16;        // 最高费率20%
    
    // 更新频率控制
    uint256 public lastUpdate;
    uint256 public constant UPDATE_INTERVAL = 1 hours;
    
    // 历史数据记录（用于分析）
    struct RateHistory {
        uint256 timestamp;
        uint256 price;
        uint256 rate;
        int256 adjustment;
    }
    
    RateHistory[] public rateHistory;
    
    function updateStabilityFee() external {
        require(block.timestamp >= lastUpdate + UPDATE_INTERVAL, "Too frequent");
        
        // 获取当前价格（可能是TWAP）
        uint256 currentPrice = priceOracle.getPrice();
        
        // 计算PID输出
        int256 adjustment = pidController.computePID(currentPrice);
        
        // 应用调整
        if (adjustment > 0) {
            currentRate = currentRate + uint256(adjustment);
        } else {
            uint256 decrease = uint256(-adjustment);
            currentRate = currentRate > decrease ? currentRate - decrease : 0;
        }
        
        // 限制范围
        if (currentRate > maxRate) {
            currentRate = maxRate;
        } else if (currentRate < minRate) {
            currentRate = minRate;
        }
        
        // 记录历史
        rateHistory.push(RateHistory({
            timestamp: block.timestamp,
            price: currentPrice,
            rate: currentRate,
            adjustment: adjustment
        }));
        
        // 更新CDP系统的费率
        cdpEngine.updateStabilityFee(currentRate);
        
        lastUpdate = block.timestamp;
        
        emit StabilityFeeUpdated(currentRate, adjustment, currentPrice);
    }
    
    // 获取系统性能指标
    function getPerformanceMetrics(uint256 periods) 
        external 
        view 
        returns (
            uint256 avgPrice,
            uint256 priceStdDev,
            uint256 avgRate,
            uint256 adjustmentCount
        ) 
    {
        require(periods > 0 && periods <= rateHistory.length, "Invalid periods");
        
        uint256 startIdx = rateHistory.length - periods;
        uint256 sumPrice;
        uint256 sumRate;
        uint256 sumSquaredDev;
        
        // 计算平均值
        for (uint256 i = startIdx; i < rateHistory.length; i++) {
            sumPrice += rateHistory[i].price;
            sumRate += rateHistory[i].rate;
            if (rateHistory[i].adjustment != 0) {
                adjustmentCount++;
            }
        }
        
        avgPrice = sumPrice / periods;
        avgRate = sumRate / periods;
        
        // 计算标准差
        for (uint256 i = startIdx; i < rateHistory.length; i++) {
            int256 deviation = int256(rateHistory[i].price) - int256(avgPrice);
            sumSquaredDev += uint256(deviation * deviation);
        }
        
        priceStdDev = sqrt(sumSquaredDev / periods);
    }
    
    // 平方根函数（用于标准差计算）
    function sqrt(uint256 x) private pure returns (uint256 y) {
        uint256 z = (x + 1) / 2;
        y = x;
        while (z < y) {
            y = z;
            z = (x / z + z) / 2;
        }
    }
}

contract LiquidationGameTheory {
    using ABDKMath64x64 for int128;
    
    struct GameParameters {
        uint256 liquidationPenalty;    // 清算罚金
        uint256 gasPrice;             // 当前Gas价格
        uint256 volatility;           // 市场波动率
        uint256 keeperCount;          // 活跃Keeper数量
    }
    
    // 计算Keeper的期望收益
    function calculateKeeperExpectedProfit(
        uint256 collateralValue,
        uint256 debtValue,
        GameParameters memory params
    ) public pure returns (int256 expectedProfit) {
        // 清算折扣收益
        uint256 discount = collateralValue * params.liquidationPenalty / 1e18;
        
        // 预期Gas成本（考虑竞争）
        uint256 gasCost = estimateGasCost(params.gasPrice, params.keeperCount);
        
        // 价格滑点风险（基于波动率）
        uint256 slippageRisk = collateralValue * params.volatility / 1e18;
        
        // 竞争风险（其他Keeper抢先）
        uint256 competitionRisk = discount * min(params.keeperCount, 10) / 20;
        
        // 期望利润 = 折扣 - Gas成本 - 滑点风险 - 竞争风险
        expectedProfit = int256(discount) - int256(gasCost) - 
                        int256(slippageRisk) - int256(competitionRisk);
    }
    
    // 计算最优抵押率（CDP持有者视角）
    function calculateOptimalCollateralRatio(
        uint256 stabilityFee,      // 年化稳定费
        uint256 collateralReturn,  // 抵押品预期收益率
        uint256 volatility,        // 波动率
        uint256 liquidationRatio   // 清算线
    ) public pure returns (uint256 optimalRatio) {
        // 使用Kelly准则的变体
        // f* = (p*b - q) / b
        // 其中：p是盈利概率，b是盈亏比，q是亏损概率
        
        // 计算持有收益
        int256 holdingReturn = int256(collateralReturn) - int256(stabilityFee);
        
        // 计算清算概率（基于波动率和距离清算线的距离）
        uint256 liquidationProb = calculateLiquidationProbability(
            liquidationRatio,
            volatility
        );
        
        // Kelly准则计算
        if (holdingReturn > 0) {
            // 有正收益时的最优抵押率
            optimalRatio = liquidationRatio * (1e18 + volatility) / 1e18;
            
            // 加入安全边际
            optimalRatio = optimalRatio * 115 / 100; // 15%安全边际
        } else {
            // 负收益时保守策略
            optimalRatio = liquidationRatio * 150 / 100;
        }
    }
    
    // 计算清算概率（使用正态分布近似）
    function calculateLiquidationProbability(
        uint256 liquidationRatio,
        uint256 annualVolatility
    ) private pure returns (uint256) {
        // 将年化波动率转换为每日波动率
        uint256 dailyVol = annualVolatility * 1e18 / sqrt(365 * 1e18);
        
        // 计算到达清算线的概率（简化模型）
        // P(清算) ≈ Φ(-d), 其中d = ln(当前比率/清算比率) / σ
        
        // 这里使用近似公式
        uint256 distance = 1e18 - (1e36 / liquidationRatio);
        uint256 probability = distance * 1e18 / (dailyVol * 3); // 简化的概率估计
        
        return probability > 1e18 ? 1e18 : probability;
    }
    
    // 协议参数优化（多目标优化）
    function optimizeProtocolParameters(
        uint256 systemRisk,
        uint256 userCost,
        uint256 keeperProfit
    ) public pure returns (
        uint256 optimalPenalty,
        uint256 optimalLiquidationRatio
    ) {
        // 目标函数：minimize(α*systemRisk + β*userCost - γ*keeperProfit)
        // 其中α、β、γ是权重参数
        
        uint256 alpha = 5;  // 系统风险权重
        uint256 beta = 3;   // 用户成本权重
        uint256 gamma = 2;  // Keeper利润权重
        
        // 使用梯度下降或其他优化算法
        // 这里简化为规则基础的调整
        
        if (systemRisk > 1e16) { // 系统风险高于1%
            optimalLiquidationRatio = 150; // 提高到150%
            optimalPenalty = 15e16;        // 15%罚金
        } else if (userCost > 5e16) { // 用户成本过高
            optimalLiquidationRatio = 130; // 降低到130%
            optimalPenalty = 10e16;        // 10%罚金
        } else {
            optimalLiquidationRatio = 140; // 标准140%
            optimalPenalty = 13e16;        // 13%罚金
        }
    }
    
    // 辅助函数
    function min(uint256 a, uint256 b) private pure returns (uint256) {
        return a < b ? a : b;
    }
    
    function sqrt(uint256 x) private pure returns (uint256) {
        if (x == 0) return 0;
        uint256 z = (x + 1) / 2;
        uint256 y = x;
        while (z < y) {
            y = z;
            z = (x / z + z) / 2;
        }
        return y;
    }
    
    function estimateGasCost(uint256 gasPrice, uint256 keeperCount) 
        private 
        pure 
        returns (uint256) 
    {
        // 基础Gas消耗
        uint256 baseGas = 200000;
        
        // 竞争导致的Gas溢价
        uint256 premiumMultiplier = 100 + keeperCount * 10;
        
        return baseGas * gasPrice * premiumMultiplier / 100;
    }
}

// 使用ABDKMath64x64库进行高精度数学计算
import "./ABDKMath64x64.sol";

contract BlackScholesCollateral {
    using ABDKMath64x64 for int128;
    
    // 常数
    int128 private constant ONE = 0x10000000000000000; // 1.0 in 64.64
    int128 private constant HALF = 0x8000000000000000; // 0.5 in 64.64
    
    struct OptionParameters {
        uint256 spotPrice;      // 抵押品现价
        uint256 strikePrice;    // 债务价值（行权价）
        uint256 timeToMaturity; // 到期时间（秒）
        uint256 volatility;     // 年化波动率（基点）
        uint256 riskFreeRate;   // 无风险利率（基点）
    }
    
    // 计算理论清算比率
    function calculateOptimalLiquidationRatio(
        OptionParameters memory params
    ) public pure returns (uint256 liquidationRatio) {
        // 转换参数为64.64定点数
        int128 S = ABDKMath64x64.divu(params.spotPrice, 1e18);
        int128 K = ABDKMath64x64.divu(params.strikePrice, 1e18);
        int128 t = ABDKMath64x64.divu(params.timeToMaturity, 365 days);
        int128 sigma = ABDKMath64x64.divu(params.volatility, 10000);
        int128 r = ABDKMath64x64.divu(params.riskFreeRate, 10000);
        
        // 计算d1和d2
        (int128 d1, int128 d2) = calculateD1D2(S, K, r, sigma, t);
        
        // 计算期权价值
        int128 putValue = calculatePutValue(S, K, r, t, d1, d2);
        
        // 清算比率 = (债务 + 期权价值) / 债务
        // 这确保了协议得到充分补偿
        int128 ratio = ONE.add(putValue.div(K));
        
        // 转换回uint256（基点）
        liquidationRatio = uint256(int256(ratio.mul(ABDKMath64x64.fromUInt(10000))));
        
        // 加入安全边际（根据波动率调整）
        uint256 safetyMargin = 10000 + params.volatility / 2; // 基础100% + 波动率的一半
        liquidationRatio = liquidationRatio * safetyMargin / 10000;
    }
    
    // 计算d1和d2
    function calculateD1D2(
        int128 S,
        int128 K,
        int128 r,
        int128 sigma,
        int128 t
    ) private pure returns (int128 d1, int128 d2) {
        // ln(S/K)
        int128 lnSK = S.div(K).ln();
        
        // (r + σ²/2)t
        int128 sigma2 = sigma.mul(sigma);
        int128 drift = r.add(sigma2.div(ABDKMath64x64.fromUInt(2))).mul(t);
        
        // σ√t
        int128 sigmaRootT = sigma.mul(t.sqrt());
        
        // d1 = [ln(S/K) + (r + σ²/2)t] / (σ√t)
        d1 = lnSK.add(drift).div(sigmaRootT);
        
        // d2 = d1 - σ√t
        d2 = d1.sub(sigmaRootT);
    }
    
    // 计算看跌期权价值
    function calculatePutValue(
        int128 S,
        int128 K,
        int128 r,
        int128 t,
        int128 d1,
        int128 d2
    ) private pure returns (int128) {
        // e^(-rt)
        int128 discountFactor = (r.neg()).mul(t).exp();
        
        // N(-d1) 和 N(-d2)
        int128 Nd1 = normalCDF(d1.neg());
        int128 Nd2 = normalCDF(d2.neg());
        
        // P = Ke^(-rt)N(-d2) - S*N(-d1)
        int128 putValue = K.mul(discountFactor).mul(Nd2).sub(S.mul(Nd1));
        
        return putValue;
    }
    
    // 正态分布累积分布函数（近似）
    function normalCDF(int128 x) private pure returns (int128) {
        // 使用Abramowitz和Stegun近似
        int128 L = x.abs();
        int128 K = ONE.div(ONE.add(L.mul(ABDKMath64x64.divu(2316419, 10000000))));
        
        int128 a1 = ABDKMath64x64.divu(319381530, 1000000000);
        int128 a2 = ABDKMath64x64.divu(356563782, 1000000000).neg();
        int128 a3 = ABDKMath64x64.divu(1781477937, 1000000000);
        int128 a4 = ABDKMath64x64.divu(1821255978, 1000000000).neg();
        int128 a5 = ABDKMath64x64.divu(1330274429, 1000000000);
        
        int128 K2 = K.mul(K);
        int128 K3 = K2.mul(K);
        int128 K4 = K3.mul(K);
        int128 K5 = K4.mul(K);
        
        int128 sum = a1.mul(K).add(a2.mul(K2)).add(a3.mul(K3))
                    .add(a4.mul(K4)).add(a5.mul(K5));
        
        // 1.0 - 1.0/√(2π) * e^(-L²/2) * sum
        int128 sqrtTwoPi = ABDKMath64x64.divu(2506628274, 1000000000);
        int128 expTerm = L.mul(L).div(ABDKMath64x64.fromUInt(2)).neg().exp();
        int128 result = ONE.sub(expTerm.div(sqrtTwoPi).mul(sum));
        
        return x >= 0 ? result : ONE.sub(result);
    }
    
    // 动态调整抵押率（基于市场条件）
    function dynamicCollateralAdjustment(
        uint256 currentVolatility,
        uint256 historicalVolatility,
        uint256 marketDepth,
        uint256 currentLiquidationRatio
    ) public pure returns (uint256 newLiquidationRatio) {
        // 波动率比率
        uint256 volRatio = currentVolatility * 1e18 / historicalVolatility;
        
        // 市场深度因子（流动性越差，需要越高的抵押率）
        uint256 depthFactor = 1e18 * 1e18 / (marketDepth + 1e18);
        
        // 调整因子 = √(波动率比率 * 深度因子)
        uint256 adjustmentFactor = sqrt(volRatio * depthFactor / 1e18);
        
        // 新抵押率 = 当前抵押率 * 调整因子
        newLiquidationRatio = currentLiquidationRatio * adjustmentFactor / 1e18;
        
        // 限制调整幅度（单次最多调整10%）
        uint256 maxChange = currentLiquidationRatio / 10;
        if (newLiquidationRatio > currentLiquidationRatio + maxChange) {
            newLiquidationRatio = currentLiquidationRatio + maxChange;
        } else if (newLiquidationRatio < currentLiquidationRatio - maxChange) {
            newLiquidationRatio = currentLiquidationRatio - maxChange;
        }
    }
    
    function sqrt(uint256 x) private pure returns (uint256) {
        if (x == 0) return 0;
        uint256 z = (x + 1) / 2;
        uint256 y = x;
        while (z < y) {
            y = z;
            z = (x / z + z) / 2;
        }
        return y;
    }
}

contract MonteCarloSimulation {
    using ABDKMath64x64 for int128;
    
    struct SimulationParameters {
        uint256 numSimulations;      // 模拟次数
        uint256 timeHorizon;        // 时间范围（天）
        uint256 initialPrice;       // 初始价格
        uint256 annualVolatility;   // 年化波动率
        uint256 drift;              // 漂移率
        uint256 jumpProbability;    // 跳跃概率
        int256 jumpSize;            // 跳跃幅度
    }
    
    struct SimulationResults {
        uint256 avgFinalPrice;      // 平均最终价格
        uint256 minPrice;           // 最低价格
        uint256 maxPrice;           // 最高价格
        uint256 defaultProbability; // 违约概率
        uint256 avgLiquidations;    // 平均清算次数
        uint256 maxDrawdown;        // 最大回撤
        uint256 var95;             // 95% VaR
        uint256 var99;             // 99% VaR
    }
    
    // 随机数生成器（简化版，实际应使用Chainlink VRF）
    uint256 private seed;
    
    function runMonteCarloSimulation(
        SimulationParameters memory params,
        uint256 liquidationRatio,
        uint256 totalCDPs
    ) public returns (SimulationResults memory results) {
        uint256[] memory finalPrices = new uint256[](params.numSimulations);
        uint256[] memory liquidationCounts = new uint256[](params.numSimulations);
        uint256 defaultCount;
        
        results.minPrice = type(uint256).max;
        
        // 运行多次模拟
        for (uint256 sim = 0; sim < params.numSimulations; sim++) {
            (uint256 finalPrice, uint256 liquidations, uint256 minPriceInPath) = 
                simulatePricePath(params, liquidationRatio);
            
            finalPrices[sim] = finalPrice;
            liquidationCounts[sim] = liquidations;
            
            // 更新统计
            results.avgFinalPrice += finalPrice;
            results.avgLiquidations += liquidations;
            
            if (finalPrice < results.minPrice) {
                results.minPrice = finalPrice;
            }
            if (finalPrice > results.maxPrice) {
                results.maxPrice = finalPrice;
            }
            
            // 计算最大回撤
            uint256 drawdown = (params.initialPrice - minPriceInPath) * 1e18 / 
                              params.initialPrice;
            if (drawdown > results.maxDrawdown) {
                results.maxDrawdown = drawdown;
            }
            
            // 检查是否违约（价格跌破某个阈值）
            if (minPriceInPath < params.initialPrice * 100 / liquidationRatio) {
                defaultCount++;
            }
        }
        
        // 计算平均值
        results.avgFinalPrice /= params.numSimulations;
        results.avgLiquidations /= params.numSimulations;
        results.defaultProbability = defaultCount * 1e18 / params.numSimulations;
        
        // 计算VaR
        sortArray(finalPrices);
        uint256 var95Index = params.numSimulations * 5 / 100;
        uint256 var99Index = params.numSimulations * 1 / 100;
        results.var95 = params.initialPrice - finalPrices[var95Index];
        results.var99 = params.initialPrice - finalPrices[var99Index];
        
        return results;
    }
    
    // 模拟单条价格路径
    function simulatePricePath(
        SimulationParameters memory params,
        uint256 liquidationRatio
    ) private returns (
        uint256 finalPrice,
        uint256 liquidationCount,
        uint256 minPrice
    ) {
        uint256 currentPrice = params.initialPrice;
        minPrice = currentPrice;
        uint256 dt = 1; // 1天时间步长
        
        // 将参数转换为每日值
        uint256 dailyVolatility = params.annualVolatility / sqrt(365);
        int256 dailyDrift = int256(params.drift) / 365;
        
        for (uint256 day = 0; day < params.timeHorizon; day++) {
            // 生成随机数
            int256 randomNormal = generateNormalRandom();
            uint256 randomUniform = generateUniformRandom();
            
            // 几何布朗运动 + 跳跃过程
            int256 diffusion = int256(dailyVolatility) * randomNormal / 1e9;
            int256 priceChange = dailyDrift + diffusion;
            
            // 检查跳跃
            if (randomUniform < params.jumpProbability) {
                priceChange += params.jumpSize;
            }
            
            // 更新价格（确保不会变负）
            if (priceChange >= 0) {
                currentPrice = currentPrice * uint256(1e18 + priceChange) / 1e18;
            } else {
                uint256 decrease = uint256(-priceChange);
                if (decrease < 1e18) {
                    currentPrice = currentPrice * (1e18 - decrease) / 1e18;
                } else {
                    currentPrice = 1; // 最小价格
                }
            }
            
            // 更新最低价
            if (currentPrice < minPrice) {
                minPrice = currentPrice;
            }
            
            // 检查清算
            if (currentPrice < params.initialPrice * 100 / liquidationRatio) {
                liquidationCount++;
            }
        }
        
        finalPrice = currentPrice;
    }
    
    // 生成正态分布随机数（Box-Muller变换）
    function generateNormalRandom() private returns (int256) {
        uint256 u1 = generateUniformRandom();
        uint256 u2 = generateUniformRandom();
        
        // Box-Muller变换
        int256 z0 = int256(sqrt(uint256(-2) * ln(u1)) * 
                    cos(2 * 314159265 * u2 / 1e10));
        
        return z0;
    }
    
    // 生成均匀分布随机数
    function generateUniformRandom() private returns (uint256) {
        seed = uint256(keccak256(abi.encodePacked(seed, block.timestamp, msg.sender)));
        return seed % 1e18;
    }
    
    // 自然对数近似
    function ln(uint256 x) private pure returns (uint256) {
        require(x > 0, "ln of non-positive number");
        // 简化的泰勒级数近似
        // 实际应用中应使用更精确的实现
        if (x == 1e18) return 0;
        
        int256 result = 0;
        int256 term = int256(x - 1e18) * 1e18 / int256(x + 1e18);
        int256 termSquared = term * term / 1e18;
        
        result = 2 * term;
        term = term * termSquared / 1e18;
        result += 2 * term / 3;
        
        return uint256(result);
    }
    
    // 余弦近似
    function cos(uint256 x) private pure returns (int256) {
        // 简化的泰勒级数
        int256 x2 = int256(x) * int256(x) / 1e18;
        int256 result = 1e18 - x2 / 2;
        return result;
    }
    
    // 数组排序（用于VaR计算）
    function sortArray(uint256[] memory arr) private pure {
        uint256 n = arr.length;
        for (uint256 i = 0; i < n - 1; i++) {
            for (uint256 j = 0; j < n - i - 1; j++) {
                if (arr[j] > arr[j + 1]) {
                    uint256 temp = arr[j];
                    arr[j] = arr[j + 1];
                    arr[j + 1] = temp;
                }
            }
        }
    }
    
    function sqrt(uint256 x) private pure returns (uint256) {
        if (x == 0) return 0;
        uint256 z = (x + 1) / 2;
        uint256 y = x;
        while (z < y) {
            y = z;
            z = (x / z + z) / 2;
        }
        return y;
    }
}

contract RiskAssessmentSystem {
    MonteCarloSimulation public simulator;
    
    struct RiskMetrics {
        uint256 systemCollateralRatio;   // 系统总抵押率
        uint256 concentrationRisk;       // 集中度风险
        uint256 liquidityRisk;          // 流动性风险
        uint256 correlationRisk;        // 相关性风险
        uint256 overallRiskScore;       // 综合风险评分
        RiskLevel riskLevel;            // 风险等级
    }
    
    enum RiskLevel { LOW, MEDIUM, HIGH, CRITICAL }
    
    // 风险阈值
    uint256 constant CRITICAL_COLLATERAL_RATIO = 120; // 120%
    uint256 constant HIGH_CONCENTRATION = 20;         // 20%单一头寸
    uint256 constant LOW_LIQUIDITY_DEPTH = 1e6 * 1e18; // 100万美元
    
    // 计算系统风险指标
    function assessSystemRisk(
        address[] memory cdps,
        address[] memory collateralTypes
    ) public returns (RiskMetrics memory metrics) {
        // 1. 计算系统总抵押率
        uint256 totalCollateralValue;
        uint256 totalDebtValue;
        
        for (uint256 i = 0; i < cdps.length; i++) {
            (uint256 collateral, uint256 debt) = getCDPValues(cdps[i]);
            totalCollateralValue += collateral;
            totalDebtValue += debt;
        }
        
        metrics.systemCollateralRatio = totalCollateralValue * 100 / totalDebtValue;
        
        // 2. 计算集中度风险（HHI指数）
        metrics.concentrationRisk = calculateHerfindahlIndex(cdps);
        
        // 3. 评估流动性风险
        metrics.liquidityRisk = assessLiquidityRisk(collateralTypes);
        
        // 4. 计算相关性风险
        metrics.correlationRisk = calculateCorrelationRisk(collateralTypes);
        
        // 5. 运行压力测试
        SimulationResults memory simResults = runStressTest();
        
        // 6. 综合评分（加权平均）
        metrics.overallRiskScore = calculateOverallRisk(
            metrics,
            simResults
        );
        
        // 7. 确定风险等级
        metrics.riskLevel = determineRiskLevel(metrics.overallRiskScore);
        
        // 8. 触发预警
        if (metrics.riskLevel >= RiskLevel.HIGH) {
            triggerRiskAlert(metrics);
        }
        
        return metrics;
    }
    
    // 计算赫芬达尔指数（市场集中度）
    function calculateHerfindahlIndex(
        address[] memory cdps
    ) private view returns (uint256) {
        uint256 totalDebt;
        uint256[] memory debts = new uint256[](cdps.length);
        
        // 获取总债务和各CDP债务
        for (uint256 i = 0; i < cdps.length; i++) {
            (, debts[i]) = getCDPValues(cdps[i]);
            totalDebt += debts[i];
        }
        
        // 计算HHI
        uint256 hhi;
        for (uint256 i = 0; i < cdps.length; i++) {
            uint256 marketShare = debts[i] * 10000 / totalDebt; // 基点
            hhi += marketShare * marketShare;
        }
        
        return hhi;
    }
    
    // 评估流动性风险
    function assessLiquidityRisk(
        address[] memory collateralTypes
    ) private view returns (uint256) {
        uint256 totalRisk;
        
        for (uint256 i = 0; i < collateralTypes.length; i++) {
            uint256 depth = getMarketDepth(collateralTypes[i]);
            
            if (depth < LOW_LIQUIDITY_DEPTH) {
                // 流动性不足，风险增加
                totalRisk += (LOW_LIQUIDITY_DEPTH - depth) * 1e18 / 
                            LOW_LIQUIDITY_DEPTH;
            }
        }
        
        return totalRisk / collateralTypes.length;
    }
    
    // 计算相关性风险
    function calculateCorrelationRisk(
        address[] memory collateralTypes
    ) private view returns (uint256) {
        // 简化版：计算价格相关系数
        // 实际应用中需要历史数据
        
        if (collateralTypes.length < 2) return 0;
        
        uint256 avgCorrelation;
        uint256 pairCount;
        
        for (uint256 i = 0; i < collateralTypes.length - 1; i++) {
            for (uint256 j = i + 1; j < collateralTypes.length; j++) {
                uint256 correlation = getPriceCorrelation(
                    collateralTypes[i],
                    collateralTypes[j]
                );
                avgCorrelation += correlation;
                pairCount++;
            }
        }
        
        return avgCorrelation / pairCount;
    }
    
    // 运行压力测试场景
    function runStressTest() private returns (SimulationResults memory) {
        SimulationParameters memory params = SimulationParameters({
            numSimulations: 1000,
            timeHorizon: 30,        // 30天
            initialPrice: 1e18,     // $1
            annualVolatility: 5000, // 50%波动率（压力场景）
            drift: -1000,          // -10%负漂移
            jumpProbability: 1e16,  // 1%跳跃概率
            jumpSize: -2e17        // -20%跳跃
        });
        
        return simulator.runMonteCarloSimulation(params, 150, 1000);
    }
    
    // 计算综合风险评分
    function calculateOverallRisk(
        RiskMetrics memory metrics,
        SimulationResults memory simResults
    ) private pure returns (uint256) {
        // 权重配置
        uint256 w1 = 30; // 抵押率权重
        uint256 w2 = 20; // 集中度权重
        uint256 w3 = 25; // 流动性权重
        uint256 w4 = 15; // 相关性权重
        uint256 w5 = 10; // 压力测试权重
        
        // 标准化各项指标（0-100）
        uint256 collateralScore = metrics.systemCollateralRatio > 200 ? 0 :
                                 (200 - metrics.systemCollateralRatio) / 2;
        
        uint256 concentrationScore = metrics.concentrationRisk / 100;
        
        uint256 liquidityScore = metrics.liquidityRisk * 100 / 1e18;
        
        uint256 correlationScore = metrics.correlationRisk * 100 / 1e18;
        
        uint256 stressScore = simResults.defaultProbability * 100 / 1e18;
        
        // 加权平均
        return (collateralScore * w1 + 
                concentrationScore * w2 + 
                liquidityScore * w3 + 
                correlationScore * w4 + 
                stressScore * w5) / 100;
    }
    
    // 确定风险等级
    function determineRiskLevel(uint256 riskScore) 
        private 
        pure 
        returns (RiskLevel) 
    {
        if (riskScore < 25) return RiskLevel.LOW;
        if (riskScore < 50) return RiskLevel.MEDIUM;
        if (riskScore < 75) return RiskLevel.HIGH;
        return RiskLevel.CRITICAL;
    }
    
    // 触发风险预警
    function triggerRiskAlert(RiskMetrics memory metrics) private {
        emit RiskAlert(
            metrics.riskLevel,
            metrics.overallRiskScore,
            block.timestamp
        );
        
        // 自动调整系统参数
        if (metrics.riskLevel == RiskLevel.CRITICAL) {
            // 提高清算线
            // 增加稳定费
            // 限制新借款
        }
    }
    
    // 辅助函数（简化版）
    function getCDPValues(address cdp) 
        private 
        view 
        returns (uint256 collateral, uint256 debt) 
    {
        // 实际实现中从CDP合约读取
        collateral = 150e18;
        debt = 100e18;
    }
    
    function getMarketDepth(address collateral) 
        private 
        view 
        returns (uint256) 
    {
        // 实际实现中从DEX获取流动性深度
        return 2e6 * 1e18;
    }
    
    function getPriceCorrelation(address token1, address token2) 
        private 
        view 
        returns (uint256) 
    {
        // 实际实现中计算历史价格相关性
        return 5e17; // 0.5相关系数
    }
    
    event RiskAlert(RiskLevel level, uint256 score, uint256 timestamp);
}