2
$\begingroup$

I am trying to simulate a Hybrid system using the WhenEvent function in Mathematica. I have asked this question earlier. It was not an MWE. Now I tried to make an MWE.

Parameter generation:

Np = 4;
classes = 2;
Xa = RandomReal[{1, 2}, {2, 2}];
Xb = RandomReal[{3, 4}, {2, 2}];
Plotxa = ListPlot[Xa, PlotStyle -> {Thick, Black}];
Plotxb = ListPlot[Xb, PlotStyle -> {Thick, Red}];
Plotx = Show[Plotxa, Plotxb, PlotRange -> {{0, 5}, {0, 5}}]
X = Join[Xa, Xb];
Xex = X /. {a_, b_} -> {1, a, b}
y = ConstantArray[1, 2]; yb = ConstantArray[-1, 2];
Y = Join[y, yb];

Below code generates the differential equations in symbolic form

β = 
  Join[{β0[t]}, Array[Symbol["β" <> ToString@#][t] &, 2]];
μ = Array[Symbol["μ" <> ToString@#][t] &, 4];
var = Join[β, μ]
varInit = var /. {t -> 0};
g = (1 - (Xex.β)*Y)
L = Total[g*μ] + 1/2 (β1[t]^2 + β2[t]^2);
f = Join[-D[L, {β}], D[L, {μ}]];

Finally the ODE simulation

DAE = {Thread[D[var, t] == f], 
   Thread[varInit == Flatten[RandomReal[{1, 2}, {1, 7}]]], 
   WhenEvent[μ1[t] == 
      0 && (1 - β0[t] - 1.0277165069495682` β1[t] - 
        1.3523023083294536` β2[t]) <= 0, 
    Derivative[1][μ1][t] -> 0]} // Flatten
sol = NDSolve[DAE, var, {t, 0, 100000}, MaxSteps -> 10000000, 
  AccuracyGoal -> 5, 
  Method -> {"EquationSimplification" -> "Residual"}]

In the DAE, I have the WhenEvent function. I was trying to model the following

$$\begin{eqnarray} \dot \mu_i&=&\left\{\begin{matrix} 0 ; &\; if\;\; \mu_i(t)=0\;\;\text{and }\;1-\beta(t)^\top x\leq0\\ 1-\beta(t)^\top x ;&\; else \end{matrix} \right. \end{eqnarray}$$

Where $x\in Xex$ set defined in the line Xex = X /. {a_, b_} -> {1, a, b}. Moreover I need to write 4 conditions as ($i\in \{1...4\}$), as a test case I was checking for just $i=1$.

From my simulation, I could see that the WhenEvent function is not working.

FYI: I added AccuracyGoal -> 5, Method -> {"EquationSimplification" -> "Residual"} because Mathematica forced me to.

$\endgroup$
5
$\begingroup$

I solved the problem without using WithEvent function but using If function. I thought I will share the solution.

Np = 4;
classes = 2;
eps = 0.01;
Xa = RandomReal[{1, 2}, {2, 2}];
Xb = RandomReal[{3, 4}, {2, 2}];
Plotxa = ListPlot[Xa, PlotStyle -> {Thick, Black}];
Plotxb = ListPlot[Xb, PlotStyle -> {Thick, Red}];
Plotx = Show[Plotxa, Plotxb, PlotRange -> {{0, 5}, {0, 5}}]
X = Join[Xa, Xb];
Xex = X /. {a_, b_} -> {1, a, b};
y = ConstantArray[1, 2]; yb = ConstantArray[-1, 2];
Y = Join[y, yb];
x = Array[Symbol["x" <> ToString@#] &, 4];
y = Array[Symbol["y" <> ToString@#] &, 4];
\[Beta] = 
  Join[{\[Beta]0[t]}, Array[Symbol["\[Beta]" <> ToString@#][t] &, 2]];
\[Mu] = Array[Symbol["\[Mu]" <> ToString@#][t] &, 4];
var = Join[\[Beta], \[Mu]]
varInit = var /. {t -> 0};
g = (1 - (Xex.\[Beta])*Y);
Cond = Thread[Thread[\[Mu] <= eps] && Thread[g <= eps]];
fmu = Transpose[{Cond, g}] /. {a_, b_} -> {If[a, 0, b]};
L = Total[g*\[Mu]] + 1/2 (\[Beta]1[t]^2 + \[Beta]2[t]^2);
f = Join[-D[L, {\[Beta]}], fmu] //  Flatten; MatrixForm[f]
opts = {Method -> {"TimeIntegration" -> {"ExplicitRungeKutta", 
       "DifferenceOrder" -> 7}}};
DAE = {Thread[D[var, t] == f], 
   Thread[varInit == Flatten[RandomReal[{100, 200}, {1, 7}]]]} // 
  Flatten; MatrixForm[DAE]
sol = NDSolve[DAE, var, {t, 0, 100000}, MaxSteps -> 10000000, 
  AccuracyGoal -> 5, 
  Method -> {"EquationSimplification" -> "Residual"}]

If someone is interested, the solution of the problem is actually an SVM (linear case).

Plot[Evaluate[(1 - (Xex.\[Beta])*Y) /. First[sol]], {t, 0, 300}, 
 PlotStyle -> {Red, Blue, Green, Orange, Black}]
Plot[Evaluate[{\[Mu]1[t], \[Mu]2[t], \[Mu]3[t], \[Mu]4[t]} /. 
   First[sol]], {t, 0, 300}, PlotStyle -> {Red, Blue, Green, Orange}]
Plot[Evaluate[{\[Beta]0[t], \[Beta]1[t], \[Beta]2[t]} /. 
   First[sol]], {t, 0, 300}, PlotStyle -> {Red, Blue, Green}]
Plotxa = ListPlot[Xa, PlotStyle -> {Thick, Black}];
Plotxb = ListPlot[Xb, PlotStyle -> {Thick, Red}];
Plotsvm = 
  Plot[-(1/\[Beta]2[t])*(\[Beta]1[t]*x + \[Beta]0[t]) /. 
     First[sol] /. {t -> 100000}, {x, 0, 2 Pi}];
Plotx = Show[Plotxa, Plotxb, Plotsvm, PlotRange -> {{0, 5}, {0, 5}}]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.