Модель САУ ЛА по лучу в среде Engee

# подключение библиотек
import Pkg
Pkg.add("ControlSystems") ## библитека для анализа САУ
Pkg.add("Plots") ## библитека для построения графиков
using ControlSystems
using Plots
using LinearAlgebra # библиотека для операций с векторами и матрицами (используется в методе градиентного спуска)
# Ограничения верхнего уровня (для расчета оптимальных коэффициентов передаточной функции)
lu = fill(300.0, 3)
# Ограничения нижнего уровня (для расчета оптимальных коэффициентов передаточной функции)
lb = fill(0.3, 3)
# Начальная точка с которой начинаем итерации по методу градиентного спуска
x0 = [0.5, 0.5, 0.5]

# Функция целей
function myfun(x::Vector{T}) where {T}
    return (150850000*x[1]) / (150850000*x[2] - 168000*x[1] + 67882500*x[3] + 49005000)
end

# Функция ограничения неравенства
function mycon(x::Vector{T}) where {T}
    c1 = 2.15 - (150850000*x[1])^-1 * (150850000*x[2] - 168000*x[1] + 67882500*x[3] + 49005000) *
              (73687500*x[3] - 150000*x[1] - 168000*x[2] + 2009250) * inv(588909 - 150000*x[2])
    
    c2 = 2.15 - (150850000*x[2] - 168000*x[1] + 67882500*x[3] + 49005000)^-1 *
              (73687500*x[3] - 150000*x[1] - 168000*x[2] + 2009250) * (588909 - 150000*x[2]) * inv(9534)
    
    c3 = 2.15 - (73687500*x[3] - 150000*x[1] - 168000*x[2] + 2009250)^-1 *
              (588909 - 150000*x[2]) * 9534 * inv(75)
    
    return [c1, c2, c3]
end

# Градиент функции цели (оцениваем численно)
function grad_f(x::Vector{T}, eps=1e-6)::Vector{T} where {T}
    n = length(x) # расчет длины вектора х
    gradient = zeros(T, n) # вектор для хранения значений градиента
    for i in 1:n
        dx = copy(x); dx[i] += eps # перебор значений х
        df = myfun(dx) - myfun(x) # расчет разности значений функции в ближайщих точках х
        gradient[i] = df / eps # расчет значений градиента
    end
    return gradient
end

# Алгоритм градиентного спуска
function minimize_gradient_descent(f, x0, max_iter=100, learning_rate=0.1, tol=1e-8)
    x = x0 # начальное значение х
    current_val = f(x0) # начальное значение функции
    prev_val = Inf  # начальное значение критерия осанова
    iter_count = 0 ## начальное значение счетчика
    while true  ## бесконечный цикл
        current_val = f(x) # текущее значение функции
        
        if abs(current_val - prev_val) < tol || iter_count >= max_iter ## критерий остановки итераций
            break
        end
        
        g = grad_f(x)  # расчет градиента функции в точке
        step = -learning_rate .* g # расчет длины шага
        
        # Проверяем границы (ограничения)
        proposed_x = clamp.(x .+ step, lb, lu) # прверка того попадает ли новая точка в интервал 
        # между заданными ограничениями
        x .= proposed_x   # обновление значения х
        prev_val = current_val  # обновление значения функции
        iter_count += 1  # увеличение счетчика итераций
    end
    return x, current_val
end

# Запускаем процесс минимизации
x_optimal, fval = minimize_gradient_descent(myfun, x0)
println("Оптимизированный вектор х:", x_optimal)
println("Значение целевой функции:", fval)

Kwz = x_optimal[1]
Kny = x_optimal[2]
Knyi = x_optimal[3]

c1 = 1.12;
c2 = 86;
c3 = 131;
c4 = 1;
c5 = 0;
c6 = 21.8;
c9 = 0.12;
g = 9.81;
a0 = c2+c1*c4;
a1 = (c1+c4+c5)
b1 = 50
b3 = 6
K = b3/b1
T = 1/b1
K1 = 1
K2 = 1

a5 = 75
a4 = 9534
a3 = (588909 - 150000*Kny)
a2 = (73687500*Kwz - 150000*Knyi - 168000*Kny + 2009250)
a1 = (150850000*Kny - 168000*Knyi + 67882500*Kwz + 49005000)
a0 = 150850000*Knyi
W2 = tf([-150000*Knyi, -168000*Knyi, 150850000*Knyi],[a5, a4, a3, a2, a1, a0])

y, t = step(W2,0:0.1:10) # расчет переходного процесса
t = collect(0:0.1:10)
plot(t,y',title="Переходный процесс") ## построение графика переходного процесса

y2, t = impulse(W2,0:0.1:10) # расчет точек импульсной характеристики
t = collect(0:0.1:10)
plot(t,y2',title="Импульсная характеристика")  # построение графика импульсной характеристики

bodeplot(W2) # построение графика частотной характеристики

nyquistplot(W2) # годограф найквиста

pzmap(W2)  # диаграмма нулей и полюсов

marginplot(W2)  # диаграмма запасов устойчивости

L = 20*log10(48)
print("Запас устойчивости: ")
print(L)
print(" Дб")