using ControlSystems
using Plots; gr()
using LinearAlgebra
A = [0 1; -4 -5]
B = [0; 1]
C = [1 0]
D = 0
P = ss(A, B, C, D)
println(P)
StateSpace{Continuous, Int64} A = 0 1 -4 -5 B = 0 1 C = 1 0 D = 0 Continuous-time state-space model
# オブザーバ極
observer_poles=[-10+5im,-10-5im]
# オブザーバゲインの設計(状態フィードバックの双対)
L = -place(P.A', P.C', observer_poles)'
println(L)
ComplexF64[-15.0 + 0.0im; -46.0 + 0.0im;;]
L = [ convert(Float64, L[i]) for i in 1:length(L) ]
2-element Vector{Float64}: -15.0 -46.0
eigvals(P.A + L * P.C)
2-element Vector{ComplexF64}: -10.0 - 5.0im -10.0 + 5.0im
G = ss(P.A, P.B, [1 0; 0 1], [0; 0]);
Obs = ss(P.A + L*P.C, [P.B -L], [1 0; 0 1], [0 0; 0 0]);
Td = 0:0.01:3;
X0 = [-1, 0.5];
# u(x,t) = [0];
Ud = zeros(1, length(Td))
y, t, x, uout = lsim(G, Ud, Td, x0=X0);
y, t, xhat, uout = lsim(Obs, [Ud' y[1,:]]', Td, [0, 0]);
p = [ plot(), plot()]
for i=1:2
plot!(p[i], t, x[i,:],
xlabel="t", #X軸のラベル
ylabel="x", #Y軸のラベル
lw=2, #線幅
ls=:solid, #線種
label="x",
size=(300,230) #プロットのサイズ
)
plot!(p[i], t, xhat[i,:],
xlabel="t", #X軸のラベル
ylabel="x", #Y軸のラベル
lw=2, #線幅
ls=:solid, #線種
label="x̂",
size=(300,230) #プロットのサイズ
)
end
plot( p[1], p[2], layout=(1,2), size=(600,230) )
# レギュレータ極
regulator_poles = [-5 + 2im, -5 - 2im]
# 極配置
F = -place(P.A, P.B, regulator_poles)
F = [convert(Float64, F[i]) for i in 1:length(F)]'
print(F)
[-25.0 -5.0]
# 出力フィードバック(オブザーバ+状態フィードバック)
K = ss(P.A+P.B*F+L*P.C, -L, F, 0)
println("K:\n", K)
println("----------------")
println("K(s)=\n", tf(K))
# フィードバック系
Gfb = feedback(P, -K);
K: StateSpace{Continuous, Float64} A = -15.0 1.0 -75.0 -10.0 B = 15.0 46.0 C = -25.0 -5.0 D = 0.0 Continuous-time state-space model ---------------- K(s)= TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}} -605.0s - 2725.0 ---------------------- 1.0s^2 + 25.0s + 225.0 Continuous-time transfer function model
Td = 0:0.01:3;
# u(x,t) = [0];
yorg, t, x, uout = lsim(P, Ud, Td, [-1, 0.5])
y, t, x, uout = lsim(Gfb, Ud, Td, [-1, 0.5, 0, 0])
plot(t, yorg',
xlabel="t", #X軸のラベル
ylabel="x", #Y軸のラベル
lw=2, #線幅
ls=:solid, #線種
legend=false,
size=(300, 230) #プロットのサイズ
)
plot!(t, y',
xlabel="t", #X軸のラベル
ylabel="x", #Y軸のラベル
lw=2, #線幅
ls=:solid, #線種
legend=false,
size=(300,230) #プロットのサイズ
)
using ControlSystems
using Plots; gr()
using LinearAlgebra
A = [0 1; -4 -5]
B = [0; 1]
C = [1 0]
D = 0
P = ss(A, B, C, D)
println(P)
A = [0 1; -4 -5]
B = [0; 1]
C = [1 0; 0 1]
D = [0; 0]
Ps = ss(A, B, C, D)
println(Ps)
StateSpace{Continuous, Int64} A = 0 1 -4 -5 B = 0 1 C = 1 0 D = 0 Continuous-time state-space model StateSpace{Continuous, Int64} A = 0 1 -4 -5 B = 0 1 C = 1 0 0 1 D = 0 0 Continuous-time state-space model
出力に0.5の定値外乱が加わるとする
# オブザーバ極
observer_poles=[-15+5im,-15-5im]
# オブザーバゲインの設計(状態フィードバックの双対)
L = -place(P.A', P.C', observer_poles)'
L = [ convert(Float64, L[i]) for i in 1:length(L) ]
G = ss(P.A , P.B, [1 0; 0 1], [0; 0]);
Obs = ss(P.A + L*P.C, [P.B -L], [1 0; 0 1], [0 0; 0 0]);
Td = 0:0.01:3;
Ud = zeros(1, length(Td))
d = zeros(1, length(Td))
d[151:length(Td)] .= 0.5
X0 = [-1, 0.5];
y, t, x, uout = lsim(G, Ud, Td, x0=X0);
y, t, xhat, uout = lsim(Obs, [Ud' y[1,:]+d']', Td, [0, 0]);
p = [ plot(), plot()]
for i=1:2
plot!(p[i], t, x[i,:],
xlabel="t", #X軸のラベル
ylabel="x", #Y軸のラベル
lw=2, #線幅
ls=:solid, #線種
label="x",
size=(300,230) #プロットのサイズ
)
plot!(p[i], t, xhat[i,:],
xlabel="t", #X軸のラベル
ylabel="x", #Y軸のラベル
lw=2, #線幅
ls=:solid, #線種
label="x̂",
size=(300,230) #プロットのサイズ
)
end
plot( p[1], p[2], layout=(1,2), size=(600,230) )
# オブザーバ極
observer_poles=[-10+5im,-10-5im, -3]
# オブザーバゲインの設計(状態フィードバックの双対)
Abar = [ P.A zeros(2,1) ; 0 0 0 ]
Bbar = [ P.B; 0 ]
Cbar = [ P.C 1 ]
Lbar = -place(Abar', Cbar', observer_poles)'
Lbar = [ convert(Float64, Lbar[i]) for i in 1:length(Lbar) ]
println(Lbar)
[75.75, -91.0, -93.75]
Aob = Abar + Lbar*Cbar
Bob = [Bbar -Lbar]
Obs = ss(Aob, Bob, [1 0 0; 0 1 0; 0 0 1], [0 0; 0 0; 0 0] );
pole(Obs)
3-element Vector{ComplexF64}: -9.999999999999982 + 5.000000000000019im -9.999999999999982 - 5.000000000000019im -3.0000000000000466 + 0.0im
Td = 0:0.01:3;
Ud = zeros(1, length(Td))
d = zeros(1, length(Td))
d[151:length(Td)] .= 0.5
X0 = [-1, 0.5];
y, t, x, uout = lsim(G, Ud, Td, x0=X0);
y, t, xhat, uout = lsim(Obs, [Ud' y[1,:]+d']', Td, [0, 0, 0]);
p = [ plot(), plot()]
for i=1:2
plot!(p[i], t, x[i,:],
xlabel="t", #X軸のラベル
ylabel="x", #Y軸のラベル
lw=2, #線幅
ls=:solid, #線種
label="x",
size=(300,230) #プロットのサイズ
)
plot!(p[i], t, xhat[i,:],
xlabel="t", #X軸のラベル
ylabel="x", #Y軸のラベル
lw=2, #線幅
ls=:solid, #線種
label="x̂",
size=(300,230) #プロットのサイズ
)
end
plot( p[1], p[2], layout=(1,2), size=(600,230) )
plot(t, xhat[3,:],
xlabel="t", #X軸のラベル
ylabel="d̂", #Y軸のラベル
lw=2, #線幅
ls=:solid, #線種
legend=false,
size=(300,230) #プロットのサイズ
)
外乱 0.5 が推定できている
using ControlSystems
using RobustAndOptimalControl
using Plots; gr()
using LinearAlgebra
g = 9.81 # 重力加速度[m/s^2]
l = 0.2 # アームの長さ[m]
M = 0.5 # アームの質量[kg]
mu = 1.5e-2 # 粘性摩擦係数
J = 1.0e-2 # 慣性モーメント
P = tf( [0,1], [J, mu, M*g*l] )
ref = 30 # 目標角度 [deg]
30
# 垂直駆動アームのノミナルモデル
Pn = tf( [0,1], [J, mu, M*g*l] )
# 不確かさ
delta = LinRange(-1, 1, 20)
WT = tf( [10, 0], [1, 150])
P = TransferFunction[ (1 + WT*delta[i])*Pn for i in 1:length(delta) ];
DT = TransferFunction[ WT*delta[i] for i in 1:length(delta) ];
setPlotScale("dB")
w = exp10.(LinRange(-2, 3, 1000))
bodeplot(Pn, w; lc=:black, lw=0.5, title="")
for sys in P
bodeplot!(sys, w; lc=:black, lw=0.5, title="" )
end
plot!(legend=false,
size=(450, 400) #プロットのサイズ
)
bodeplot(DT[1], w; lc=:black, lw=0.5, title="")
for sys in DT
bodeplot!(sys, w; lc=:black, lw=0.5, title="" )
end
plot!(legend=false,
size=(450, 400) #プロットのサイズ
)
確認中
WS = tf( [0, 1], [1, 1, 0.25]); # 感度関数に対する重み関数
WU = ss(1.0)
WT = tf( [10, 0], [1, 150]); # 相補感度関数に対する重み関数
G = hinfpartition(Pn, WS, WU, WT)
# Check if the system is feasible for synythesis
flag = hinfassumptions(G, verbose=false)
true
# Synthesize the H-infinity optimal controller
K, γ = hinfsynthesize(G, γrel=1.0);
Pcl, S, CS, T = hinfsignals(G, Pn, K);
println(tf(K))
println(γ)
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}} 6.914688563198183s^4 + 1053.488413717986s^3 + 3129.9699783801843s^2 + 103660.1636521904s + 87011.21342993349 ----------------------------------------------------------------------------------------------------------------------- 1.0s^5 + 164.802459568436s^4 + 2400.3155738379014s^3 + 24135.781867094975s^2 + 22417.633916500767s + 5464.6420344833505 Continuous-time transfer function model 0.9527801586107133
w = exp10.(LinRange(-3, 3, 1000))
# 感度関数
Ssys = feedback(1, Pn * K)
bodeplot(Ssys, w; plotphase=false, lw=2, title="")
bodeplot!(γ / WS, w; plotphase=false, lc=:black, title="")
# 相補感度関数
Tsys = feedback(Pn * K, 1)
bodeplot!(Tsys, w; plotphase=false, lw=2, title="")
bodeplot!(γ / WT, w; plotphase=false, lc=:black, title="")
plot!(legend=false,
ylims=(-60, 40),
size=(450, 200) #プロットのサイズ
)
bodeplot(K, w; plotphase=false, lw=2, title="")
bodeplot!(Pn, w; plotphase=false, lw=2, title="")
plot!(legend=false,
ylims=(-80, 40),
size=(450, 200) #プロットのサイズ
)
## Plot the specifications
specificationplot([S, CS, T], [WS, WU, WT], γ) |> display
plt = plot()
for i in 1:1:length(delta)
P = (1 + WT * delta[i]) * Pn
Gyr = feedback(P * K, 1)
y, t = step(Gyr, 0:0.01:2)
plot!(plt, t, y', lc=:black, lw=0.3)
end
plot!(plt,
xlabel="t", #X軸のラベル
ylabel="y", #Y軸のラベル
ylim=(0, 2),
lw=2, #線幅
ls=:solid, #線種
legend=false,
size=(300, 230) #プロットのサイズ
)
H = TransferFunction[K * (1 + WT * delta[i]) * Pn for i in 1:length(delta)];
nyquistplot(H[1], linecolor=:black, linewidth=1)
for sys in H
nyquistplot!(sys, linecolor=:black, linewidth=1, legend=false)
end
plot!(xlims=(-1.5, 0.5),
ylims=(-1.5, 0.5),
size=(300, 300), #プロットのサイズ
)
┌ Warning: Keyword argument hover not supported with Plots.GRBackend(). Choose from: annotationcolor, annotationfontfamily, annotationfontsize, annotationhalign, annotationrotation, annotations, annotationvalign, arrow, aspect_ratio, axis, background_color, background_color_inside, background_color_outside, background_color_subplot, bar_width, bins, bottom_margin, camera, clims, color_palette, colorbar, colorbar_entry, colorbar_scale, colorbar_title, colorbar_titlefont, colorbar_titlefontcolor, colorbar_titlefontrotation, colorbar_titlefontsize, connections, contour_labels, discrete_values, fill, fill_z, fillalpha, fillcolor, fillrange, fillstyle, flip, fontfamily, fontfamily_subplot, foreground_color, foreground_color_axis, foreground_color_border, foreground_color_grid, foreground_color_subplot, foreground_color_text, formatter, framestyle, grid, gridalpha, gridlinewidth, gridstyle, group, guide, guidefont, guidefontcolor, guidefontfamily, guidefonthalign, guidefontrotation, guidefontsize, guidefontvalign, html_output_format, inset_subplots, label, layout, left_margin, legend_background_color, legend_column, legend_font, legend_font_color, legend_font_family, legend_font_halign, legend_font_pointsize, legend_font_rotation, legend_font_valign, legend_foreground_color, legend_position, legend_title, legend_title_font_color, legend_title_font_family, legend_title_font_pointsize, legend_title_font_rotation, legend_title_font_valigm, levels, lims, line, line_z, linealpha, linecolor, linestyle, linewidth, link, margin, marker_z, markeralpha, markercolor, markershape, markersize, markerstrokealpha, markerstrokecolor, markerstrokewidth, minorgrid, minorgridalpha, minorgridlinewidth, minorgridstyle, minorticks, mirror, normalize, orientation, overwrite_figure, permute, plot_title, plot_titlefontcolor, plot_titlefontfamily, plot_titlefontrotation, plot_titlefontsize, plot_titlelocation, plot_titlevspan, polar, primary, projection, quiver, ribbon, right_margin, rotation, scale, series_annotations, seriesalpha, seriescolor, seriestype, show, show_empty_bins, showaxis, size, smooth, subplot, subplot_index, thickness_scaling, tick_direction, tickfontcolor, tickfontfamily, tickfonthalign, tickfontrotation, tickfontsize, tickfontvalign, ticks, title, titlefontcolor, titlefontfamily, titlefonthalign, titlefontrotation, titlefontsize, titlefontvalign, top_margin, unitformat, weights, widen, window_title, x, xdiscrete_values, xerror, xflip, xforeground_color_axis, xforeground_color_border, xforeground_color_grid, xforeground_color_text, xformatter, xgrid, xgridalpha, xgridlinewidth, xgridstyle, xguide, xguidefontcolor, xguidefontfamily, xguidefonthalign, xguidefontrotation, xguidefontsize, xguidefontvalign, xlims, xlink, xminorgrid, xminorgridalpha, xminorgridlinewidth, xminorgridstyle, xminorticks, xmirror, xrotation, xscale, xshowaxis, xtick_direction, xtickfontcolor, xtickfontfamily, xtickfonthalign, xtickfontrotation, xtickfontsize, xtickfontvalign, xticks, xunitformat, xwiden, y, ydiscrete_values, yerror, yflip, yforeground_color_axis, yforeground_color_border, yforeground_color_grid, yforeground_color_text, yformatter, ygrid, ygridalpha, ygridlinewidth, ygridstyle, yguide, yguidefontcolor, yguidefontfamily, yguidefonthalign, yguidefontrotation, yguidefontsize, yguidefontvalign, ylims, ylink, yminorgrid, yminorgridalpha, yminorgridlinewidth, yminorgridstyle, yminorticks, ymirror, yrotation, yscale, yshowaxis, ytick_direction, ytickfontcolor, ytickfontfamily, ytickfonthalign, ytickfontrotation, ytickfontsize, ytickfontvalign, yticks, yunitformat, ywiden, z, z_order, zdiscrete_values, zerror, zflip, zforeground_color_axis, zforeground_color_border, zforeground_color_grid, zforeground_color_text, zformatter, zgrid, zgridalpha, zgridlinewidth, zgridstyle, zguide, zguidefontcolor, zguidefontfamily, zguidefonthalign, zguidefontrotation, zguidefontsize, zguidefontvalign, zlims, zlink, zminorgrid, zminorgridalpha, zminorgridlinewidth, zminorgridstyle, zminorticks, zmirror, zrotation, zscale, zshowaxis, ztick_direction, ztickfontcolor, ztickfontfamily, ztickfonthalign, ztickfontrotation, ztickfontsize, ztickfontvalign, zticks, zunitformat, zwiden └ @ Plots ~/.julia/packages/Plots/sxUvK/src/args.jl:1557
kp = 2
kd = 0.1
ki = 10
plt = plot()
for i in 1:1:length(delta)
Kpid = tf([kd, kp, ki], [1, 0])
P = (1 + WT*delta[i])*Pn
Gyr = feedback(P*Kpid, 1)
y, t = step(Gyr, 0:0.01:2 )
plot!(plt, t, y', lc=:black, lw =0.3)
end
plot!(plt,
xlabel="t", #X軸のラベル
ylabel="y", #Y軸のラベル
ylim = (0,2),
lw=2, #線幅
ls=:solid, #線種
legend=false,
size=(300,230) #プロットのサイズ
)
ナイキスト線図で不確かさの影響を確認する
Kpid = tf([kd, kp, ki], [1, 0])
H = TransferFunction[ Kpid*(1 + WT*delta[i])*Pn for i in 1:length(delta) ];
nyquistplot(H[1], linecolor=:black, linewidth=1)
for sys in H
nyquistplot!(sys, linecolor=:black, linewidth=1, legend=false)
end
plot!(xlims=(-1.5, 0.5),
ylims=(-1.5, 0.5),
size=(300, 300), #プロットのサイズ
)
using ControlSystems
using Plots; gr()
using LinearAlgebra
P = tf([0, 1], [0.5, 1])
println("連続時間システム",P)
連続時間システムTransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}} 1.0 ---------- 0.5s + 1.0 Continuous-time transfer function model
ts = 0.2
Pd1 = c2d(P, ts, :zoh)
println("離散時間システム(zoh)", Pd1)
離散時間システム(zoh)TransferFunction{Discrete{Float64}, ControlSystemsBase.SisoRational{Float64}} 0.3296799539643608 ------------------------- 1.0z - 0.6703200460356393 Sample Time: 0.2 (seconds) Discrete-time transfer function model
Pd2 = c2d(P, ts, :tustin)
print("離散時間システム(tustin)",Pd2)
離散時間システム(tustin)TransferFunction{Discrete{Float64}, ControlSystemsBase.SisoRational{Float64}} 0.16666666666666669z + 0.1666666666666667 ----------------------------------------- 1.0z - 0.6666666666666667 Sample Time: 0.2 (seconds) Discrete-time transfer function model
Tc = 0:0.01:3
y, t = step(P, Tc)
Td = 0:ts:3;
p = [ plot(), plot()]
plot!(p[1], t, y',
xlabel="t", #X軸のラベル
ylabel="y", #Y軸のラベル
lw=2, #線幅
ls=:solid, #線種
label="continuous",
size=(300,230) #プロットのサイズ
)
y, t = step(Pd1, Td)
plot!(p[1], t, y',
seriestype = :scatter,
xlabel="t", #X軸のラベル
ylabel="y", #Y軸のラベル
label="zoh",
legend=:bottomright,
size=(300,230) #プロットのサイズ
)
plot!(p[2], t, y',
xlabel="t", #X軸のラベル
ylabel="y", #Y軸のラベル
lw=2, #線幅
ls=:solid, #線種
label="continuous",
size=(300,230) #プロットのサイズ
)
y, t = step(Pd2, Td)
plot!(p[2], t, y',
seriestype = :scatter,
xlabel="t", #X軸のラベル
ylabel="y", #Y軸のラベル
label="tustin",
legend=:bottomright,
size=(300,230) #プロットのサイズ
)
plot( p[1], p[2], layout=(1,2), size=(600,230) )
Tc = 0:0.01:3
u(x,t) = [0.5 * sin.(6*t) + 0.5 * cos.(8*t)]
yc, tc, x, uout = lsim(P, u, Tc);
Td = 0:ts:3;
p = [ plot(), plot()]
plot!(p[1], tc, yc',
xlabel="t", #X軸のラベル
ylabel="y", #Y軸のラベル
lw=2, #線幅
ls=:solid, #線種
label="continuous",
size=(300,230) #プロットのサイズ
)
y, t, x, uout = lsim(Pd1, u, Td);
plot!(p[1], t, y',
seriestype = :scatter,
xlabel="t", #X軸のラベル
ylabel="y", #Y軸のラベル
label="zoh",
legend=:bottomleft,
size=(300,230) #プロットのサイズ
)
plot!(p[2], tc, yc',
xlabel="t", #X軸のラベル
ylabel="y", #Y軸のラベル
lw=2, #線幅
ls=:solid, #線種
label="continuous",
size=(300,230) #プロットのサイズ
)
y, t, x, uout = lsim(Pd2, u, Td);
plot!(p[2], t, y',
seriestype = :scatter,
xlabel="t", #X軸のラベル
ylabel="y", #Y軸のラベル
label="tustin",
legend=:bottomleft,
size=(300,230) #プロットのサイズ
)
plot( p[1], p[2], layout=(1,2), size=(600,230) )
setPlotScale("dB")
bodeplot(P, lw=2, size=(450,400), legend=false, title="" )
setPlotScale("dB")
bodeplot([Pd1, Pd2], lw=2, size=(450,400), legend=false, title="" )
using ControlSystems
using JuMP
using Ipopt # 非線形最適化ソルバー
using Plots
# パラメータ
xf = 0.0 # 基準状態
uf = 0.0 # 基準入力
Q = 50.0 # 状態に関する重み
R = 0.1 # 入力に関する重み
umin, umax = -4.0, 4.0 # 入力の制約
Tf = 0.5 # 最終時間
dt = 0.01 # 時間ステップ
N = Int(Tf / dt) # ステップ数
# 最適化モデル
model = Model(Ipopt.Optimizer)
@variables(model, begin
-10 <= x[0:N] <= 10 # 状態変数
umin <= uopt[0:N-1] <= umax # 制御入力
end)
# 目的関数(二次コスト)
@objective(model, Min, sum(Q * (x[i] - xf)^2 + R * (uopt[i] - uf)^2 for i in 0:N-1))
# 動的制約(状態方程式)
@constraint(model, [i in 0:N-2], x[i+1] == x[i] + dt * (x[i] + uopt[i]))
# 初期状態
@constraint(model, x[0] == 1.0)
# 最適化問題を解く
optimize!(model)
# 結果の出力
println("最適化成功: ", termination_status(model))
# println("最適な状態: ", value.(x))
# println("最適な入力: ", value.(u))
****************************************************************************** This program contains Ipopt, a library for large-scale nonlinear optimization. Ipopt is released as open source code under the Eclipse Public License (EPL). For more information visit https://github.com/coin-or/Ipopt ****************************************************************************** This is Ipopt version 3.14.13, running with linear solver MUMPS 5.6.2. Number of nonzeros in equality constraint Jacobian...: 148 Number of nonzeros in inequality constraint Jacobian.: 0 Number of nonzeros in Lagrangian Hessian.............: 100 Total number of variables............................: 101 variables with only lower bounds: 0 variables with lower and upper bounds: 101 variables with only upper bounds: 0 Total number of equality constraints.................: 50 Total number of inequality constraints...............: 0 inequality constraints with only lower bounds: 0 inequality constraints with lower and upper bounds: 0 inequality constraints with only upper bounds: 0 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 0 0.0000000e+00 1.00e+00 0.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0 1 3.1967994e+01 6.77e-01 4.78e-01 -1.0 1.23e+01 - 2.45e-01 3.23e-01h 1 2 6.6563619e+01 5.53e-01 2.29e+00 -1.0 5.70e+00 - 5.23e-01 1.84e-01h 1 3 1.7523577e+02 3.44e-01 2.92e+00 -1.0 4.51e+00 - 6.15e-01 3.77e-01h 1 4 3.0964589e+02 1.95e-01 3.75e+00 -1.0 3.44e+00 - 8.71e-01 4.32e-01h 1 5 4.4437684e+02 8.93e-02 2.54e+00 -1.0 2.48e+00 - 1.00e+00 5.43e-01h 1 6 5.4551095e+02 2.47e-02 8.31e-01 -1.0 1.46e+00 - 1.00e+00 7.23e-01h 1 7 5.8680596e+02 1.53e-16 1.81e-13 -1.7 4.93e-01 - 1.00e+00 1.00e+00h 1 8 5.8639232e+02 1.73e-16 2.11e-13 -2.5 1.32e-01 - 1.00e+00 1.00e+00f 1 9 5.8631681e+02 1.04e-16 1.09e-13 -3.8 9.34e-02 - 1.00e+00 1.00e+00f 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 10 5.8631520e+02 1.25e-16 1.69e-13 -3.8 4.97e-02 - 1.00e+00 1.00e+00f 1 11 5.8631137e+02 1.73e-16 1.50e-13 -5.7 2.72e-02 - 1.00e+00 1.00e+00f 1 12 5.8631130e+02 1.67e-16 1.37e-13 -5.7 1.36e-02 - 1.00e+00 1.00e+00f 1 13 5.8631129e+02 1.67e-16 2.63e-13 -5.7 6.49e-03 - 1.00e+00 1.00e+00f 1 14 5.8631124e+02 1.46e-16 3.38e-13 -8.6 3.15e-03 - 1.00e+00 1.00e+00f 1 15 5.8631124e+02 1.04e-16 2.00e-13 -8.6 1.07e-03 - 1.00e+00 1.00e+00f 1 16 5.8631124e+02 6.25e-17 3.16e-13 -8.6 1.59e-04 - 1.00e+00 1.00e+00h 1 Number of Iterations....: 16 (scaled) (unscaled) Objective...............: 5.8631124411209680e+02 5.8631124411209680e+02 Dual infeasibility......: 3.1572071873383807e-13 3.1572071873383807e-13 Constraint violation....: 6.2450045135165055e-17 6.2450045135165055e-17 Variable bound violation: 3.9841570043108732e-08 3.9841570043108732e-08 Complementarity.........: 8.8875886325376534e-09 8.8875886325376534e-09 Overall NLP error.......: 8.8875886325376534e-09 8.8875886325376534e-09 Number of objective function evaluations = 17 Number of objective gradient evaluations = 17 Number of equality constraint evaluations = 17 Number of inequality constraint evaluations = 0 Number of equality constraint Jacobian evaluations = 1 Number of inequality constraint Jacobian evaluations = 0 Number of Lagrangian Hessian evaluations = 1 Total seconds in IPOPT = 0.015 EXIT: Optimal Solution Found. 最適化成功: LOCALLY_SOLVED
# システムの定義
A, B, C, D = 1, 1, 1, 0
sys = ss(A, B, C, D)
# シミュレーション時間と初期状態
Td = 0:0.01:0.5
x0 = [1.0]
# 最適化結果から入力を取得
u_optimal = [value(uopt[i]) for i in 0:N-1]
push!(u_optimal, u_optimal[end])
u_optimal_vec = vec(u_optimal)
# システムの応答を計算
y, t, x, uout = lsim(sys, u_optimal_vec', Td, x0)
# 終端状態を表示
println("Final simulated state: ", x[:, end])
# プロット
p1 = plot(t, x', label="x", lw=2, xlabel="t", ylabel="x")
p2 = plot(t, uout', label="u", lw=2, xlabel="t", ylabel="u")
plot(p1, p2, layout=(1, 2), size=(600, 230), legend=false)
Final simulated state: [-0.0056711182693998195]