(* by Nasser M. Abbasi, version: April 16, 2012*)
Manipulate[
tick;
needToTick = False;
currentTime = currentTime + delT;
If[( (currentTime + delT) <= simulationTime) && (run || oneStep),
(
invertedPendulum@makeStep[delT];
If[oneStep,
(
oneStep = False;
run = False
),
(
needToTick = True
)
]
)];
{cx, ctheta, cxSpeed, cThetaSpeed} =
invertedPendulum@getCurrentPosition[];
display@addPoint[cx, cxSpeed, ctheta, cThetaSpeed, currentTime];
finalPlot = display@getPlot[];
If[needToTick, tick += del];
FinishDynamic[];
finalPlot,
Grid[{
{Button[
Text@Style["run", 12], {run = True; oneStep = False;
display@reset[simulationTime, delT, ic]; currentTime = 0;
tick += del}, ImageSize -> {50, 35}],
Button[
Text@Style["step", 12], {run = True; oneStep = True,
tick += del}, ImageSize -> {50, 35}]},
{
Button[
Text@Style["stop", 12], {run = False; oneStep = False;
tick += del}, ImageSize -> {50, 35}],
Button[
Text@Style["reset", 12], {currentTime = -delT;
invertedPendulum@reset[];
display@reset[simulationTime, delT, ic]; run = False;
oneStep = False; tick += del}, ImageSize -> {50, 35}]
}
}, Spacings -> {0.2, 0.2}, Alignment -> Center],
Grid[{
{Text@Style["time", 12],
Manipulator[
Dynamic[simulationTime, {simulationTime = #; currentTime = -delT;
invertedPendulum@reset[];
display@reset[simulationTime, delT, ic]; run = False;
tick += del} &], {1, 100, 1}, ImageSize -> Tiny],
Text@Style[Dynamic@padIt2[simulationTime, {3, 0}], 11]
},
{Text@Style["slow", 12],
Manipulator[
Dynamic[delT, {delT = #; currentTime = -delT;
invertedPendulum@reset[];
display@reset[simulationTime, delT, ic];
run = False} &], {0.01, 0.3, 0.01}, ImageSize -> Tiny],
Text@Style["fast", 12]
}
}, Frame -> True, FrameStyle -> Directive[Thickness[.001], Gray]
],
Grid[{
{Text@Style["\[Theta](0)", 12],
Manipulator[
Dynamic[ic, {ic = #; currentTime = -delT; run = False;
display@reset[simulationTime, delT, ic];
invertedPendulum@setInitialAngle[ic*Pi/180.];
tick += del} &], {70, 110, 1}, ImageSize -> Tiny],
Text[Row[{Style[Dynamic@padIt2[ic, {3, 0}], 11], Degree}]]
},
{Text@Style[Row[{Style["x", Italic], "(0)"}], 12],
Manipulator[
Dynamic[icx, {icx = #; currentTime = -delT; run = False;
display@reset[simulationTime, delT, ic];
invertedPendulum@setInitialX[icx]; tick += del} &], {-2, 2,
0.1}, ImageSize -> Tiny],
Text[Row[{Style[Dynamic@padIt1[icx, {2, 1}], 11]}]]
}
}, Frame -> True, FrameStyle -> Directive[Thickness[.001], Gray]],
Grid[{
{Text@Style["bob mass", 12],
Manipulator[
Dynamic[bobMass, {bobMass = #; currentTime = -delT; run = False;
display@reset[simulationTime, delT, ic];
invertedPendulum@setBobMass[bobMass]; tick += del} &], {0.1,
10, 0.1}, ImageSize -> Tiny],
Text@Style[Dynamic@padIt2[bobMass, {3, 1}], 11]
},
{Text@Style["cart mass", 12],
Manipulator[
Dynamic[cartMass, {cartMass = #; currentTime = -delT;
run = False; display@reset[simulationTime, delT, ic];
invertedPendulum@setCartMass[cartMass]; tick += del} &], {0.1,
10, 0.1}, ImageSize -> Tiny],
Text@Style[Dynamic@padIt2[cartMass, {3, 1}], 11]
},
{Text@Style["length", 12],
Manipulator[
Dynamic[pendulumLength, {pendulumLength = #;
currentTime = -delT; run = False;
display@reset[simulationTime, delT, ic];
invertedPendulum@setPendulumLength[pendulumLength];
tick += del} &], {1, 2, 0.01}, ImageSize -> Tiny],
Text@Style[Dynamic@padIt2[pendulumLength, {2, 1}], 11]
}
}, Frame -> True, FrameStyle -> Directive[Thickness[.001], Gray]
],
Grid[{
{Text@Style["Q matrix diagonal", 12], SpanFromLeft},
{Text@Style[
Row[{Style["x", Italic], "(", Style["t", Italic], ")"}], 12],
Manipulator[
Dynamic[positionWeight, {positionWeight = #;
currentTime = -delT; run = False;
display@reset[simulationTime, delT, ic];
invertedPendulum@setPositionWeight[positionWeight];
tick += del} &], {1, 1000, 1}, ImageSize -> Tiny],
Text@Style[Dynamic@padIt2[positionWeight, {4, 0}], 11]
},
{Text@Style[
Row[{Style["x'", Italic], "(", Style["t", Italic], ")"}], 12],
Manipulator[
Dynamic[linearVelocityWeight, {linearVelocityWeight = #;
currentTime = -delT; run = False;
display@reset[simulationTime, delT, ic];
invertedPendulum@
setLinearVelocityWeight[linearVelocityWeight];
tick += del} &], {1, 1000, 1}, ImageSize -> Tiny],
Text@Style[Dynamic@padIt2[linearVelocityWeight, {4, 0}], 12]
},
{Text@Style[Row[{"\[Theta](", Style["t", Italic], ")"}], 12],
Manipulator[
Dynamic[angleWeight, {angleWeight = #; currentTime = -delT;
run = False; display@reset[simulationTime, delT, ic];
invertedPendulum@setAngleWeight[angleWeight];
tick += del} &], {1, 1000, 1}, ImageSize -> Tiny],
Text@Style[Dynamic@padIt2[angleWeight, {4, 0}], 11]
},
{Text@Style[
Row[{"\[Theta]\[InvisibleSpace]'(", Style["t", Italic], ")"}],
12],
Manipulator[
Dynamic[angluarVelocityWeight, {angluarVelocityWeight = #;
currentTime = -delT; run = False;
display@reset[simulationTime, delT, ic];
invertedPendulum@
setAngluarVelocityWeight[angluarVelocityWeight];
tick += del} &], {1, 1000, 1}, ImageSize -> Tiny],
Text@Style[Dynamic@padIt2[angluarVelocityWeight, {4, 0}], 11]
}}, Frame -> True, FrameStyle -> Directive[Thickness[.001], Gray]
],
Grid[{
{Text@Style["friction coefficients", 12], SpanFromLeft},
{Text@Style["static", 11],
Manipulator[
Dynamic[staticFrictionCoefficient, {staticFrictionCoefficient = \
#; currentTime = -delT; run = False;
invertedPendulum@setStaticFriction[staticFrictionCoefficient];
display@reset[simulationTime, delT, ic]; tick += del} &], {0,
0.1, 0.05}, ImageSize -> Tiny],
Text@Style[Dynamic@padIt2[staticFrictionCoefficient, {3, 2}], 11]
},
{Text@Style["kinetic", 11],
Manipulator[
Dynamic[kineticFrictionCoefficient, {kineticFrictionCoefficient \
= #; currentTime = -delT; run = False;
invertedPendulum@
setKineticFriction[kineticFrictionCoefficient];
display@reset[simulationTime, delT, ic]; tick += del} &], {0,
0.05, 0.01}, ImageSize -> Tiny],
Text@Style[Dynamic@padIt2[kineticFrictionCoefficient, {3, 2}], 11]
},
{Text@Style["viscous", 11],
Manipulator[
Dynamic[viscousFrictionCoefficient, {viscousFrictionCoefficient \
= #; currentTime = -delT; run = False;
invertedPendulum@
setViscousFriction[viscousFrictionCoefficient];
display@reset[simulationTime, delT, ic]; tick += del} &], {0,
6, 0.05}, ImageSize -> Tiny],
Text@Style[Dynamic@padIt2[viscousFrictionCoefficient, {3, 2}], 11]
}
}, Alignment -> Center, Frame -> True,
FrameStyle -> Directive[Thickness[.001], Gray]
],
{{oneStep, False}, None},
{{currentTime, -0.1}, None},
{tick, None},
{{run, False}, None},
{p, None},
{cx, None},
{ctheta, None},
{cxSpeed, None},
{cThetaSpeed, None},
{{ic, 75}, None},
{{icx, 1}, None},
{{bobMass, 1}, None},
{{cartMass, 10}, None},
{{pendulumLength, 2}, None},
{{delT, 0.1}, None},
{{del, $MachineEpsilon}, None},
{finalPlot, None},
{{positionWeight, 100}, None},
{{linearVelocityWeight, 10}, None},
{{angleWeight, 10}, None},
{{angluarVelocityWeight, 1}, None},
{{simulationTime, 10}, None},
{{staticFrictionCoefficient, 0.05}, None},
{{kineticFrictionCoefficient, 0.03}, None},
{{viscousFrictionCoefficient, 3.0}, None},
{needToTick, None},
TrackedSymbols :> {tick},
SynchronousUpdating -> False,
ContinuousAction -> False,
SynchronousInitialization -> True,
Alignment -> Center,
ImageMargins -> 0,
FrameMargins -> 0,
Paneled -> True,
Frame -> False,
ControlPlacement -> Left,
AutorunSequencing -> {1},
Initialization :>
{
ContentSizeW = 430;
ContentSizeH = 520 ;
(*-----------------------------------------*)
padIt1[v_?(NumericQ[#] && Im[#] == 0 &), f_List] :=
AccountingForm[Chop[N@v] , f, NumberSigns -> {"-", "+"},
NumberPadding -> {"0", "0"}, SignPadding -> True];
(*-----------------------------------------*)
padIt2[v_?(NumericQ[#] && Im[#] == 0 &), f_List] :=
AccountingForm[Chop[N@v] , f, NumberSigns -> {"", ""},
NumberPadding -> {"0", "0"}, SignPadding -> True];
(*-----------------------------------------*)
invertedPendulumClass[$bobMass_, $cartMass_, $pendulumLen_, \
$staticFriction_, $kineticFriction_, $viscousFriction_, \
$initialAngle_, $initialX_, $stepSize_, $positionWeight_,
$linearVelocityWeight_, $angleWeight_, $angluarVelocityWeight_, \
$x_, $\[Theta]_, $t_] := Module[{
bobMass = $bobMass,
cartMass = $cartMass,
pendulumLen = $pendulumLen,
staticFriction = $staticFriction,
kineticFriction = $kineticFriction,
viscousFriction = $viscousFriction,
initialAngle = $initialAngle,
initialX = $initialX,
positionWeight = $positionWeight,
linearVelocityWeight = $linearVelocityWeight,
angleWeight = $angleWeight,
angluarVelocityWeight = $angluarVelocityWeight,
x = $x,
\[Theta] = $\[Theta],
t = $t,
eqns,
lastAngle = $initialAngle,
lastX = $initialX,
lastSpeed = 0,
lastAngularSpeed = 0,
lastAppliedForce = 0,
lastCoulombForce = 0,
lastViscousForce = 0,
normalForce,
stateControlExpression,
closedLoopEigenvalues,
gain,
f,
model,
self,
update,
init,
calculateFrictionForce},
(*private methods *)
(*----------------------------------------*)
init[] :=
Module[{ke, pe, lag, eq1, eq2, g = 9.8, currentSolution},
ke =
1/2 cartMass x'[t]^2 +
1/2 bobMass * (x'[t]^2 + pendulumLen^2*\[Theta]'[t]^2 -
2 x'[t]*pendulumLen*\[Theta]'[t] Sin[\[Theta][t]]);
pe = bobMass * g * pendulumLen* Sin[\[Theta][t]];
lag = ke - pe;
eq1 = D[D[lag, \[Theta]'[t]], t] - D[lag, \[Theta][t]] == 0;
eq2 =
D[D[lag, x'[t]], t] - D[lag, x[t]] ==
f[t] - viscousFriction*x'[t];
eqns = {eq1, eq2};
model =
StateSpaceModel[
eqns, {{x[t], 0}, {x'[t], 0}, {\[Theta][t],
Pi/2}, {\[Theta]'[t], 0}}, {{f[t], 0}}, {\[Theta][t],
x[t]}, t];
gain =
First@LQRegulatorGains[
N[model], {DiagonalMatrix[{positionWeight,
linearVelocityWeight, angleWeight,
angluarVelocityWeight}], {{1}}}];
stateControlExpression = {x[t],
x'[t], \[Theta][t] - Pi/2, \[Theta]'[t]};
closedLoopEigenvalues = {{Eigenvalues[
First[Normal[
SystemsModelStateFeedbackConnect[N[model], {gain}]]]]}};
lastAngle = initialAngle;
lastX = initialX;
lastSpeed = 0;
lastAngularSpeed = 0;
normalForce = (cartMass + bobMass)*9.8
];
(*public methods*)
(*----------------------------------------*)
self@makeStep[delT_] :=
Module[{currentSolution, initialConditions, effectiveForce},
initialConditions = {x[0] == lastX,
x'[0] == lastSpeed, \[Theta]'[0] ==
lastAngularSpeed, \[Theta][0] == lastAngle};
lastAppliedForce = -gain. stateControlExpression /. {x[t] ->
lastX, \[Theta][t] -> lastAngle,
x'[t] -> lastSpeed, \[Theta]'[t] -> lastAngularSpeed};
currentSolution =
First@If[Abs[lastSpeed] <= $MachineEpsilon,
(
effectiveForce =
If[Abs[lastAppliedForce] < normalForce*staticFriction,
0,
-(gain.stateControlExpression) -
normalForce*staticFriction*Sign[lastAppliedForce]
];
NDSolve[
Join[eqns /. f[t] -> (effectiveForce*UnitStep[t]),
initialConditions], {x, \[Theta], x', \[Theta]'}, {t, 0,
delT}]
),
(
lastCoulombForce = -kineticFriction*normalForce*
Sign[lastSpeed];
lastViscousForce = -viscousFriction*lastSpeed;
NDSolve[
Join[eqns /.
f[t] -> (-(gain.stateControlExpression) +
lastCoulombForce*UnitStep[t]),
initialConditions], {x, \[Theta], x', \[Theta]'}, {t, 0,
delT}]
)
];
lastSpeed = (x' /. currentSolution)[delT];
lastAngle = (\[Theta] /. currentSolution)[delT];
lastX = (x /. currentSolution)[delT];
lastAngularSpeed = (\[Theta]' /. currentSolution)[delT]
];
self@reset[] := (init[]);
self@setStaticFriction[v_] := (staticFriction = v; init[]);
self@setKineticFriction[v_] := (kineticFriction = v; init[]);
self@setViscousFriction[v_] := (viscousFriction = v; init[]);
self@setCartMass[v_] := (cartMass = v; init[]);
self@setBobMass[v_] := (bobMass = v; init[]);
self@setPendulumLength[v_] := (pendulumLen = v; init[]);
self@setInitialAngle[v_] := (initialAngle = v; init[]);
self@setInitialX[v_] := (initialX = v; init[]);
self@setPositionWeight[v_] := (positionWeight = v; init[]);
self@setLinearVelocityWeight[v_] := (linearVelocityWeight = v;
init[]);
self@setAngleWeight[v_] := (angleWeight = v; init[]);
self@setAngluarVelocityWeight[v_] := (angluarVelocityWeight = v;
init[]);
(*----------------------------------------*)
self@
getCurrentPosition[] := ({lastX, lastAngle, lastSpeed,
lastAngularSpeed});
self@getGain[] := gain;
self@getFrictionType[] := (
Text@
Style[If[Abs[lastSpeed] <= $MachineEpsilon, "static",
"kinetic"], 11]);
self@getStateSpace[] := model;
self@getAppliedForce[] := lastAppliedForce;
self@getCoulombForce[] := lastCoulombForce;
self@getViscousForce[] := lastViscousForce;
self@getPendulumLength[] := pendulumLen;
self@getClosedLoopEigenvalues[] := closedLoopEigenvalues;
(*constructor*)
init[];
self
];
(*----------------------------------------*)
displayClass[$simulationTime_, $delT_, $initialAngle_] :=
Module[{simulationTime = $simulationTime, delT = $delT,
initialAngle = $initialAngle, xyData, currentIndex, init,
makePolesAndZerosCoordinates, polesPlot, getStatistics, self},
(*----------------------------------------*)
init[] := (
xyData =
Table[{0, 0, 0, 0, 0}, {Floor[simulationTime/delT] + 1}];
currentIndex = 0;
xyData[[1, 4]] = initialAngle*Pi/180.
);
(*-----------------------------------------*)
makePolesAndZerosCoordinates[points_] := Module[{xy},
xy =
Flatten[Replace[
points, {Complex[x_, y_] :> {x, y},
x_?NumericQ :> {x, 0}}, {3}], 1];
Cases[xy, {_?NumericQ, _?NumericQ}, {2}]
];
(*-----------------------------------------*)
polesPlot[poles_, plotStyle_] :=
Module[{polesPoints, plotOptions},
plotOptions = {AxesOrigin -> {0, 0},
ImagePadding -> {{20, 20}, {20, 0}},
Frame -> True,
ImageMargins -> 1, AspectRatio -> 1};
polesPoints = makePolesAndZerosCoordinates[poles];
ListPlot[polesPoints,
AxesOrigin -> {0, 0},
PlotMarkers -> Style["\[CircleTimes]", plotStyle, 11],
GridLines -> Automatic,
GridLinesStyle ->
Directive[Dashed, Thickness[.001], LightGray],
Evaluate@plotOptions,
PlotRange -> {
{ Min[-1.2, 1.2*Min[polesPoints[[All, 1]] ] ],
Max[1.2, 1.2*Max[polesPoints[[All, 1]]] ]}, {
Min[-1.2, 1.2*Min[polesPoints[[All, 2]] ] ],
Max[1.2, 1.2*Max[polesPoints[[All, 2]]] ]}
},
ImageSize -> {0.4 ContentSizeW, 0.4 ContentSizeH},
ImageMargins -> 0,
PlotLabel ->
Text[Row[{"[", Style["A-BK", Italic, 12],
Style["] eigenvalues", 12]}]]
]
];
(*----------------------------------------*)
self@reset[$simulationTime2_, $delT2_, $initialAngle2_] :=
(
simulationTime = $simulationTime2;
delT = $delT2;
initialAngle = $initialAngle2;
init[]
);
(*----------------------------------------*)
self@
addPoint[x_, xSpeed_, angle_, angleSpeed_, currentTime_] :=
(
currentIndex++;
xyData[[currentIndex]] = {currentTime, x, xSpeed, angle*180/Pi,
angleSpeed}
);
(*----------------------------------------*)
getStatistics[] := Module[{c},
If[currentIndex == 0, c = 1, c = currentIndex];
Grid[{
{Text@
Style[Row[{"time = ", padIt2[xyData[[c, 1]], {4, 2}],
" sec"}], 12],
Text[Row[{Style["x", Italic, 12], " = ",
Style[padIt1[xyData[[c, 2]], {4, 3}], 12], " ",
Style["m", Italic, 12]}]],
Text@Style[
Row[{"\[Theta] = ", padIt2[xyData[[c, 4]], {4, 2}],
" deg"}], 12]
},
{
Text[Row[{Style["x'(t)", Italic, 12], " = ",
padIt1[xyData[[c, 3]], {5, 4}], Style[" m/sec", 12]}]],
Text[Style[
Row[{"\[Theta] '(t) = ", padIt1[xyData[[c, 5]], {5, 4}],
" rad/sec"}], 12]],
""
}
,
{Text[
Row[{Style["K", Italic, 12],
Style[" (state gain vector) = ", 12],
TraditionalForm@
NumberForm[
Style[invertedPendulum@getGain[], 11], {5, 2},
SignPadding -> True, NumberPadding -> {"0", "0"},
NumberSigns -> {"-", "+"}]}]], SpanFromLeft
},
{Text[
Row[{Style["f", Italic, 12], " = ",
Style[padIt1[invertedPendulum@getAppliedForce[], {5, 3}],
12], " ", Style["N", Italic, 12]}]],
Text[Row[{Style["\!\(\*SubscriptBox[\(F\), \(c\)]\) ",
Italic, 12], " = ",
Style[padIt1[invertedPendulum@getCoulombForce[], {5, 3}],
12], " ", Style["N", Italic, 12]}]],
Text[Row[{Style["\!\(\*SubscriptBox[\(F\), \(v\)]\)",
Italic, 12], " = ",
Style[padIt1[invertedPendulum@getViscousForce[], {5, 3}],
12], " ", Style["N", Italic, 12]}]]
}
}, Frame -> All, Spacings -> {0.5, .2},
FrameStyle -> Directive[Thickness[.001], Gray]
]
];
(*----------------------------------------*)
self@getPlot[] := Module[{g0, g1, g2, g3, g4, c, len},
g0 =
polesPlot[invertedPendulum@getClosedLoopEigenvalues[], Blue];
len = invertedPendulum@getPendulumLength[];
If[currentIndex == 0, c = 1, c = currentIndex];
g1 = Graphics[
{
{Blue,
Rectangle[{xyData[[c, 2]] - 1/2, 0}, {xyData[[c, 2]] + 1/2,
3/15}]},
{
Red,
Line[{
{xyData[[c, 2]], 3/15},
{xyData[[c, 2]] + len*Cos[xyData[[c, 4]]*Pi/180],
3/15 + len*Sin[xyData[[c, 4]]*Pi/180]}
}]
},
{Disk[{xyData[[c, 2]] + len*Cos[xyData[[c, 4]]*Pi/180],
3/15 + len*Sin[xyData[[c, 4]]*Pi/180]}, .15]}
},
PlotRange -> {{-4, 4}, {-0.5, 1.4*len}},
ImageSize -> {0.4 ContentSizeW, 0.4 ContentSizeH},
AspectRatio -> 1,
Frame -> True,
Axes -> True,
AxesStyle -> Dashed,
FrameLabel -> {{None, None}, {Style["x", Italic, 12], None}},
ImagePadding -> {{22, 5}, {33, 5}},
ImageMargins -> 0
];
g2 = ListPlot[xyData[[1 ;; c, 1 ;; 2]],
PlotRange -> {{0, simulationTime}, All},
Joined -> True,
ImagePadding -> {{40, 10}, {36, 35}},
Frame -> True,
ImageSize -> {0.51 ContentSizeW, 0.43 ContentSizeH},
FrameLabel -> {{None, None}, {Text@Style["time (sec)", 12],
Text@Style[" cart position vs. time", 12]}},
ImageMargins -> 0,
GridLines -> Automatic,
GridLinesStyle ->
Directive[Dashed, Thickness[.001], LightGray],
PlotStyle -> Red,
AspectRatio -> 1,
Axes -> None,
Epilog -> {Dashed, Thin, Line[{{0, 0}, {simulationTime, 0}}]}
];
g3 = ListPlot[xyData[[1 ;; c, {1, 4}]],
PlotRange -> {{0, simulationTime}, All},
Joined -> True,
ImagePadding -> {{40, 10}, {36, 35}},
Frame -> True,
ImageSize -> {0.51 ContentSizeW, 0.43 ContentSizeH},
FrameLabel -> {{None, None}, {Text@Style["time (sec)", 12],
Text@Style[" bob angle vs. time", 12]}},
ImageMargins -> 0,
GridLines -> Automatic,
GridLinesStyle ->
Directive[Dashed, Thickness[.001], LightGray],
PlotStyle -> Red,
AspectRatio -> 1,
Axes -> None,
Epilog -> {Dashed, Thin,
Line[{{0, 90}, {simulationTime, 90}}]}
];
g4 = getStatistics[];
Grid[{{g4, SpanFromLeft}, {g0, g1}, {g2, g3}},
Spacings -> {0, 0}, Alignment -> Center, Frame -> All,
FrameStyle -> Directive[Thickness[.001], Gray]]
];
init[];
self
];
invertedPendulum =
invertedPendulumClass[1, 10, 2, 0.1, 0.03, 3, 75*Pi/180., 1, 0.05,
100, 10, 10, 1, x, \[Theta], t];
display = displayClass[10, 0.05, 5];
}
]