已知系统的开环传递函数为 G(s)H(s)=k(s+1)/s(s-1)(s^2+4s+20)试确定使系统阻尼比为0.5的k值和这时的闭环特征根。 （注意：不能用目测定位。） 写出matlab简单易懂代码

在控制系统理论中，给定开环传递函数 \( G(s) \cdot H(s) = \frac{k(s+1)}{s(s-1)(s^2 + 4s + 20)} \)，我们想要找到使得系统阻尼比 \( \zeta \) 等于 0.5 的 \( k \) 值，并计算此时的闭环极点。

对于一阶二型系统 (Type I system)，当阻尼比 \( \zeta \) 等于 0.5 时，系统是一个临界稳定的系统，即两个实部相等的复数极点对称分布在虚轴上。因此，我们需要解方程来找出 \( k \) 的值，使得 \( s = -\omega_n \pm j\omega_n \)，其中 \( \omega_n \) 是自然频率，\( j \) 是虚数单位。

对于给定的 \( G(s)H(s) \)，我们可以尝试用MATLAB来求解。首先，我们需要找到 \( s \) 满足条件的表达式，然后代入到传递函数中，找到对应的 \( k \) 值。

% 定义系数
a1 = 1;
b1 = 1;
c1 = 0;
d1 = 0;
a2 = 0;
b2 = 1;
c2 = 4;
d2 = 20;

% 自然频率 ωn 的近似值，因为无法直接目测，这里假设一个初始值
wn_guess = sqrt(c2 / d2); % 自然频率

% 极点的初始估计，考虑阻尼比为0.5的情况
zn_guess = 0.5;
p = wn_guess * [1; zn_guess];

% 使用 MATLAB 的 roots 函数查找闭环极点
K = 1; % 初始设定 k=1 作为搜索起点
while true
    K_guess = K;
    % 闭合环传递函数
    GHS = @(s) K_guess * a1*s + b1*(s+1) ./ (s.*(s-1).*((s+p(1)).*(s+p(2))));
    
    % 求解闭环传递函数的零点
    poles_closed_loop = roots(GHS);
    
    % 计算阻尼比
    zeta = real(poles_closed_loop(1)) / wn_guess;
    
    % 如果阻尼比接近目标值，停止循环
    if abs(zeta - 0.5) < eps
        break;
    end
    
    % 更新 k 值并继续搜索
    K = K_guess / (real(poles_closed_loop(1))^2);
end

% 输出结果
fprintf('为了达到阻尼比为 0.5，k 值应为: %.6f\n', K);
fprintf('闭环特征根为: %.6f + j%.6f 和 %.6f - j%.6f\n', poles_closed_loop(1), imag(poles_closed_loop(1)), poles_closed_loop(2), -imag(poles_closed_loop(2)));

这个代码会不断迭代更新 \( K \) 值直到阻尼比达到 0.5。请注意，由于阻尼比的计算依赖于数值解法，实际运行时可能会有误差，尤其是在初始猜测 \( \omega_n \) 或 \( K \) 不精确的情况下。如果需要更精确的结果，可以使用数值优化算法如 fminsearch。