I implement an LQR control for a three-story structure.

Firstly, I achieved a reasonable result with my first code (`Code 1`

As shown below). In this code, the output mainly contains the shifts of three floors.

In my second code (`Code 2`

If I use a different output that contains mostly ground accelerations, Mathematica returns the error "General" :: "munfl". Thanks to the help of Suba Thomas, I've found that state variables can change if I want an output acceleration. I think I used a proper state variable in my second code.

Another confusing thing is that my second code returns a different result (see screenshot below), even though I used the same output as in the first code.

My questions are:

- Why my second code gives completely different results, though

Was the same output used?
- How can I successfully solve Mathematica in my second code by using an output that differs from the output in Code 1?

The data file "elcentro_NS.dat" used in my codes can be downloaded from http://www.vibrationdata.com/elcentro.htm.

Suggestions are welcome.

Thank you very much 🙂

`Code 1:`

(gives reasonable results)

```
Remove("Global`*") // Quiet;
SetDirectory(NotebookDirectory());
SetOptions(Plot, PlotRange -> All, ImageSize -> Large,
Evaluated -> True, Frame -> True );
txgDat = Cases(
Import("elcentro_NS.dat", "Data"), {x1_, x2_} :> {x1, 9.8 x2});
tDat = txgDat((All, 1));
tend = Max(tDat);
xg = Interpolation(txgDat)(t);
xgMax = Max(Abs@txgDat((All, 2)));
pNmax = xgMax/1000;
mNmax = xgMax/1000;
w = 0.1;
v = 0.01;
pDat = {tDat,
pNmax RandomVariate(NormalDistribution(0, Sqrt(w)),
Length(tDat))}(Transpose);
mDat = {tDat,
mNmax RandomVariate(NormalDistribution(0, Sqrt(v)),
Length(tDat))}(Transpose);
pNoise = Interpolation(pDat)(t);
mNoise = Interpolation(mDat)(t);
(* -- m, k symbolic matrix for finding correct state ------ *)
(*m={{Subscript((ScriptM), 1),0,0},{0,Subscript((ScriptM),
2),0},{0,0,Subscript((ScriptM), 3)}};
k=Array(Subscript((ScriptK), #1,#2)&,{3,3});
c=Array(Subscript((ScriptC), #1,#2)&,{3,3});*)
(*
----------------------------------------------------- *)
(* -------------- m, k numerical matrix ----------------- *)
m = {{m1, 0, 0}, {0, m2, 0}, {0, 0, m3}};
k = {{k1 + k2, -k2, 0}, {-k2, k2 + k3, -k3}, {0, -k3, k3}};
m1 = m2 = m3 = 4 10^5;
k1 = k2 = k3 = 2 10^8;
(Zeta)1 = (Zeta)2 = 5/100;
(Omega)sol = (Omega) /.
NSolve({Det(k - (Omega)^2 m) == 0, (Omega) > 0}, (Omega));
{(Omega)1, (Omega)2} = (Omega)sol((1 ;; 2));
(Alpha)c = (
2 (Omega)1 (Omega)2 ((Zeta)1 (Omega)2 - (Zeta)2 (Omega)1))/(
(Omega)2^2 - (Omega)1^2);
(Beta)c = (
2 ((Zeta)2 (Omega)2 - (Zeta)1 (Omega)1))/((Omega)2^2 -
(Omega)1^2);
c = (Alpha)c m + (Beta)c k; c/10^6;
(* ---------------------------------------------------- *)
bs = {{1, -1, 0}, {0, 1, -1}, {0, 0, 1}};
(CapitalLambda) = {{1, 1, 1}}(Transpose);
uf = {uf1(t), uf2(t), uf3(t)};(*反馈输入*)
ue = {ue1(t), ue2(t), ue3(t)};
pN = {pN1(t), pN2(t), pN3(t)};
mN = {mN1(t), mN2(t), mN3(t)};
u = Flatten@{uf, ue, pN, mN};
eq = Flatten(
m.{x1''(t), x2''(t), x3''(t)} + c.{x1'(t), x2'(t), x3'(t)} +
k.{x1(t), x2(t), x3(t)}) ==
Flatten(-m.ue + bs.uf - m.(CapitalLambda) pN) // Thread;
(DoubleStruckY) = {x1(t), x2(t), x3(t)} + ue + mN;
(DoubleStruckZ) = {x1(t), x2(t), x3(t), x1'(t), x2'(t), x3'(t)};
ss = StateSpaceModel(eq, (DoubleStruckZ), u, (DoubleStruckY), t,
SystemsModelLabels -> {ToString /@ u, ToString /@ y,
ToString /@ z});
zUncontrol =
StateResponse(
ss, {0, 0, 0, xg, xg, xg, pNoise, pNoise, pNoise, mNoise, mNoise,
mNoise}, {t, 0, tend})((1 ;; 3));
(*Grid@{Table(Plot(zUncontrol(LeftDoubleBracket)i
(RightDoubleBracket),{t,0,tend},PlotTheme(Rule)"Web"),{i,3})}*)
(Alpha) = 50; (Beta) = 5 10^-6;
qMatrix = (Alpha) ArrayFlatten({{k, 0}, {0,
m}})(*IdentityMatrix(6)*);
rMatrix = (Beta) IdentityMatrix(3);
kLQR = LQRegulatorGains({ss, {1, 2, 3}}, {qMatrix, rMatrix});
Flqr = -kLQR.(DoubleStruckZ);
zLQRode = NDSolveValue(
{eq /.
{uf1(t) -> Flqr((1)), uf2(t) -> Flqr((2)),
uf3(t) -> Flqr((3)),
ue1(t) -> xg, ue2(t) -> xg, ue3(t) -> xg,
pN1(t) -> pNoise, pN2(t) -> pNoise, pN3(t) -> pNoise,
mN1(t) -> mNoise, mN2(t) -> mNoise, mN3(t) -> mNoise} //
Simplify,
x1(0) == x2(0) == x3(0) == x1'(0) == x2'(0) == x3'(0) == 0}, {x1(
t), x2(t), x3(t)}, {t, 0, tend});
(*Grid@{Table(Plot(zLQRode(LeftDoubleBracket)i(RightDoubleBracket),{
t,0,30},PlotTheme(Rule)"Web"),{i,3})}*)
Table(Plot({zUncontrol((i)), zLQRode((i)), PlotTheme -> "Web"}, {t, 0,
tend}), {i, 3})
```

`Code 2:`

(returns no or false result if the same result as in code 1 is used)

```
Remove("Global`*") // Quiet;
SetDirectory(NotebookDirectory());
SetOptions(Plot, PlotRange -> All, ImageSize -> Large,
Evaluated -> True, Frame -> True );
txgDat = Cases(
Import("elcentro_NS.dat", "Data"), {x1_, x2_} :> {x1, 9.8 x2});
tDat = txgDat((All, 1));
tend = Max(tDat);
xg = Interpolation(txgDat)(t);
xgMax = Max(Abs@txgDat((All, 2)));
pNmax = xgMax/1000;
mNmax = xgMax/1000;
w = 0.1;
v = 0.01;
pDat = {tDat,
pNmax RandomVariate(NormalDistribution(0, Sqrt(w)),
Length(tDat))}(Transpose);
mDat = {tDat,
mNmax RandomVariate(NormalDistribution(0, Sqrt(v)),
Length(tDat))}(Transpose);
pNoise = Interpolation(pDat)(t);
mNoise = Interpolation(mDat)(t);
Plot({pNoise, mNoise}, {t, 0, tend}, PlotRange -> All,
Evaluated -> True);
(* -- m, k symbolic matrix for finding correct state ------ *)
(*m={{Subscript((ScriptM), 1),0,0},{0,Subscript((ScriptM),
2),0},{0,0,Subscript((ScriptM), 3)}};
k=Array(Subscript((ScriptK), #1,#2)&,{3,3});
c=Array(Subscript((ScriptC), #1,#2)&,{3,3});*)
(*
----------------------------------------------------- *)
(* -------------- m, k numerical matrix ----------------- *)
m = {{m1, 0, 0}, {0, m2, 0}, {0, 0, m3}};
k = {{k1 + k2, -k2, 0}, {-k2, k2 + k3, -k3}, {0, -k3, k3}};
m1 = m2 = m3 = 4 10^5;
k1 = k2 = k3 = 2 10^8;
(Zeta)1 = (Zeta)2 = 5/100;
(Omega)sol = (Omega) /.
NSolve({Det(k - (Omega)^2 m) == 0, (Omega) > 0}, (Omega));
{(Omega)1, (Omega)2} = (Omega)sol((1 ;; 2));
(Alpha)c = (
2 (Omega)1 (Omega)2 ((Zeta)1 (Omega)2 - (Zeta)2 (Omega)1))/(
(Omega)2^2 - (Omega)1^2);
(Beta)c = (
2 ((Zeta)2 (Omega)2 - (Zeta)1 (Omega)1))/((Omega)2^2 -
(Omega)1^2);
c = (Alpha)c m + (Beta)c k; c/10^6;
(* ---------------------------------------------------- *)
bs = {{1, -1, 0}, {0, 1, -1}, {0, 0, 1}};
(CapitalLambda) = {{1, 1, 1}}(Transpose);
uf = {uf1(t), uf2(t), uf3(t)};
ue = {ue1(t), ue2(t), ue3(t)};
pN = {pN1(t), pN2(t), pN3(t)};
mN = {mN1(t), mN2(t), mN3(t)};
u = Flatten@{uf, ue, pN, mN};
eq = Flatten(
m.{x1''(t), x2''(t), x3''(t)} + c.{x1'(t), x2'(t), x3'(t)} +
k.{x1(t), x2(t), x3(t)}) ==
Flatten(-m.ue + bs.uf - m.(CapitalLambda) pN) // Thread;
y = {x1''(t), x2''(t), x3''(t)} + ue +
mN;(* with this y, error of "General::munfl" generates *)
(*y={x1(t),x2(t),x3(t)}+ue+mN;*) (* with this y, wrong
result generates *)
z = {x3'(t), x3(t), x2'(t), x2(t), x1'(t), x1(t)};
(* state changed from descriptor ssm to a standard ssm *)
ss = StateSpaceModel(eq, z, u, y, t,
SystemsModelLabels -> {ToString /@ u, ToString /@ y,
ToString /@ z});
zUncontrol =
StateResponse(
ss, {0, 0, 0, xg, xg, xg, pNoise, pNoise, pNoise, mNoise, mNoise,
mNoise}, {t, 0, tend})(({6, 4, 2}));
(*Grid@{Table(Plot(zUncontrol(LeftDoubleBracket)i
(RightDoubleBracket),{t,0,tend},PlotTheme(Rule)"Web"),{i,3})}*)
(Alpha) = 50; (Beta) = 5 10^-6;
qMatrix = (Alpha) ArrayFlatten({{k, 0}, {0,
m}})(*IdentityMatrix(6)*);
rMatrix = (Beta) IdentityMatrix(3);
kLQR = LQRegulatorGains({ss, {1, 2, 3}}, {qMatrix, rMatrix});
Flqr = -kLQR.z;
zLQRode = NDSolveValue(
{eq /.
{uf1(t) -> Flqr((1)), uf2(t) -> Flqr((2)),
uf3(t) -> Flqr((3)),
ue1(t) -> xg, ue2(t) -> xg, ue3(t) -> xg,
pN1(t) -> pNoise, pN2(t) -> pNoise, pN3(t) -> pNoise,
mN1(t) -> mNoise, mN2(t) -> mNoise, mN3(t) -> mNoise} //
Simplify,
x1(0) == x2(0) == x3(0) == x1'(0) == x2'(0) == x3'(0) == 0}, {x1(
t), x2(t), x3(t)}, {t, 0, tend});
(*Grid@{Table(Plot(zLQRode(LeftDoubleBracket)i(RightDoubleBracket),{
t,0,30},PlotTheme(Rule)"Web"),{i,3})}*)
Table(Plot({zUncontrol((i)), zLQRode((i))}, {t, 0, tend}), {i, 3})
```