tf, ss mldivide

nelson-lang · Oct 17, 2023 · 6f52a2e · 6f52a2e
1 parent 04fca70
commit 6f52a2e
Show file tree

Hide file tree

Showing 6 changed files with 174 additions and 50 deletions.
diff --git a/modules/control_system/functions/@ss/mldivide.m b/modules/control_system/functions/@ss/mldivide.m
@@ -0,0 +1,18 @@
+%=============================================================================
+% Copyright (c) 2023-present Allan CORNET (Nelson)
+%=============================================================================
+% This file is part of the Nelson.
+%=============================================================================
+% LICENCE_BLOCK_BEGIN
+% SPDX-License-Identifier: LGPL-3.0-or-later
+% LICENCE_BLOCK_END
+%=============================================================================
+function varargout = mldivide(varargin)
+    narginchk(2, 2),
+    nargoutchk(0, 1);
+    sysA = ss(varargin{1});
+    sysB = ss(varargin{2});
+    sys = inv(sysA) * sysB;
+    varargout{1} = sys;
+end
+%=============================================================================
diff --git a/modules/control_system/functions/@ss/mtimes.m b/modules/control_system/functions/@ss/mtimes.m
@@ -8,63 +8,80 @@
 % LICENCE_BLOCK_END
 %=============================================================================
 function sys = mtimes(sys1, sys2)
-    sysA = ss(sys1);
-    sysB = ss(sys2);
-
-    Ts = mtimes_timesample(sysA.Ts, sysB.Ts);
-    if isempty(sysA.A) && isempty(sysA.B) && isempty(sysA.C) && isscalar(sysA.D)
-        A = sysB.A;
-        B = sysA.D * sysB.B; 
-        C = sysB.C;
-        D = sysA.D * sysB.D;
+  sysA = ss(sys1);
+  sysB = ss(sys2);
+
+  Ts = mtimes_timesample(sysA.Ts, sysB.Ts);
+  if isempty(sysA.A) && isempty(sysA.B) && isempty(sysA.C) && isscalar(sysA.D)
+    A = sysB.A;
+    B = sysA.D * sysB.B; 
+    C = sysB.C;
+    D = sysA.D * sysB.D;
+  else
+    if isempty(sysB.A) && isempty(sysB.B) && isempty(sysB.C) && isscalar(sysB.D)
+      A = sysA.A;
+      B = sysB.D * sysA.B;
+      C = sysA.C;
+      D = sysA.D * sysB.D;
     else
-        if isempty(sysB.A) && isempty(sysB.B) && isempty(sysB.C) && isscalar(sysB.D)
-            A = sysA.A;
-            B = sysB.D * sysA.B;
-            C = sysA.C;
-            D = sysA.D * sysB.D;
-        else
-            [pA, mA] = size(sysA);
-            [pB, mB] = size(sysB);
-            if (mA ~= pA)
-                error(_('Matrix dimensions must agree.'));
-            end
-            A = [[sysA.A, sysA.B * sysB.C];[zeros(size(sysB.A, 1), size(sysA.A,2)), sysB.A]];
-            B = [sysA.B*sysB.D;sysB.B];
-            C = [sysA.C, sysA.D * sysB.C];
-            D = sysA.D * sysB.D;
-        end
-    end
-    sys = ss(A, B, C, D, Ts);
-    UserData = [];
-    if ~isempty(sysA.UserData) && ~isempty(sysB.UserData)
-      UserData = sysB.UserData;
-    else
-      if ~isempty(sysA.UserData)
-        UserData = sysA.UserData;
-      end
-      if ~isempty(sysA.UserData)
-        UserData = sysA.UserData;
+      [pA, mA] = size(sysA);
+      [pB, mB] = size(sysB);
+      if (mA ~= pA)
+        error(_('Matrix dimensions must agree.'));
       end
+      A = [[sysA.A, sysA.B * sysB.C];[zeros(size(sysB.A, 1), size(sysA.A,2)), sysB.A]];
+      B = [sysA.B*sysB.D;sysB.B];
+      C = [sysA.C, sysA.D * sysB.C];
+      D = sysA.D * sysB.D;
     end
-    if ~isempty(UserData)
-      sys.UserData = UserData;
+  end
+
+  sys = ss(A, B, C, D, Ts);
+  sys.E = computesE(sysA.E,sysB.E, size(sysA.A, 1),size(sysB.A,1));
+  UserData = [];
+  if ~isempty(sysA.UserData) && ~isempty(sysB.UserData)
+    UserData = sysB.UserData;
+  else
+    if ~isempty(sysA.UserData)
+      UserData = sysA.UserData;
     end
+    if ~isempty(sysA.UserData)
+      UserData = sysA.UserData;
+    end
+  end
+  if ~isempty(UserData)
+    sys.UserData = UserData;
+  end
 end
 %=============================================================================
 function Ts = mtimes_timesample(TsA, TsB)
-    if (TsA == -1) && (TsB ~= -1)
-        Ts = TsB;
+  if (TsA == -1) && (TsB ~= -1)
+    Ts = TsB;
+  else
+    if (TsB == -1) && (TsA ~= -1)
+      Ts = TsA;
     else
-        if (TsB == -1) && (TsA ~= -1)
-            Ts = TsA;
-        else
-            if (TsB == TsA)
-                Ts = TsA;
-            else
-                error(_('Sampling times must agree.'));
-            end
-        end
-    end 
+      if (TsB == TsA)
+        Ts = TsA;
+      else
+        error(_('Sampling times must agree.'));
+      end
+    end
+  end 
+end
+%=============================================================================
+function E = computesE(E1, E2, n1, n2)
+  if ~isempty(E1) || ~isempty(E2)
+    if isempty(E1)
+      E1 = eye(n1,'like', E2);
+    elseif isempty(E2)
+      E2 = eye(n2,'like', E1);
+    end
+    E = zeros(n1+n2, n1+n2);
+    E(1:n1,1:n1) = E1;
+    E(n1+1:n1 + n2,n1+1:n1+n2) = E2; 
+  else
+    E = [];
+  end
 end
 %=============================================================================
diff --git a/modules/control_system/functions/@tf/mldivide.m b/modules/control_system/functions/@tf/mldivide.m
@@ -0,0 +1,18 @@
+%=============================================================================
+% Copyright (c) 2023-present Allan CORNET (Nelson)
+%=============================================================================
+% This file is part of the Nelson.
+%=============================================================================
+% LICENCE_BLOCK_BEGIN
+% SPDX-License-Identifier: LGPL-3.0-or-later
+% LICENCE_BLOCK_END
+%=============================================================================
+function varargout = mldivide(varargin)
+    narginchk(2, 2),
+    nargoutchk(0, 1);
+    sysA = tf(varargin{1});
+    sysB = tf(varargin{2});
+    sys = inv(sysA) * sysB;
+    varargout{1} = sys;
+end
+%=============================================================================
diff --git a/modules/control_system/tests/test_ss_mldivide.m b/modules/control_system/tests/test_ss_mldivide.m
@@ -0,0 +1,43 @@
+%=============================================================================
+% Copyright (c) 2023-present Allan CORNET (Nelson)
+%=============================================================================
+% This file is part of the Nelson.
+%=============================================================================
+% LICENCE_BLOCK_BEGIN
+% SPDX-License-Identifier: LGPL-3.0-or-later
+% LICENCE_BLOCK_END
+%=============================================================================
+A = [1 2; 3 4];
+B = [1 0; 0 1];
+C = [1 1; 1 1];
+D = [1 2; 3 4];
+sys1 = ss(A, B, C, D);
+sys2 = ss(A, B, C, D);
+R = sys1 \ sys2;
+
+REF_A = [1   2   1   0   0   0;
+ 3   4   0   1   0   0;
+ 1   1   1   2  -1  -1;
+ 1   1   3   4  -1  -1;
+ 0   0   0   0   1   2;
+ 0   0   0   0   3   4];
+REF_B = [0   0;
+0   0;
+-1  -2;
+-3  -4;
+1   0;
+0   1];
+REF_C = [0   0   1   0   0   0;
+0   0   0   1   0   0];
+REF_D = [0   0;
+0   0];
+REF_E = [1   0   0   0   0   0;
+         0   1   0   0   0   0;
+         0   0   0   0   0   0;
+         0   0   0   0   0   0;
+         0   0   0   0   1   0;
+         0   0   0   0   0   1];
+REF = ss(REF_A, REF_B, REF_C, REF_D);
+REF.E = REF_E;
+assert_isequal(R, REF);
+
diff --git a/modules/control_system/tests/test_tf_mldivide.m b/modules/control_system/tests/test_tf_mldivide.m
@@ -0,0 +1,25 @@
+%=============================================================================
+% Copyright (c) 2023-present Allan CORNET (Nelson)
+%=============================================================================
+% This file is part of the Nelson.
+%=============================================================================
+% LICENCE_BLOCK_BEGIN
+% SPDX-License-Identifier: LGPL-3.0-or-later
+% LICENCE_BLOCK_END
+%=============================================================================
+sysA = tf(3);
+sysB = tf(4);
+sys = sysA \ sysB;
+REF = tf(4/3);
+assert_isequal(sys, REF)
+%=============================================================================
+num = [3 4];
+den = [3 1 5];
+Ts = 0.2;
+sys = tf(num, den, Ts);
+R = sys \ sys;
+num_REF = [9  15  19  20];
+den_REF = [9  15  19  20];
+REF = tf(num_REF, den_REF, Ts);
+assert_isequal(R, REF);
+%=============================================================================
diff --git a/modules/elementary_functions/src/cpp/IsEqualTo.cpp b/modules/elementary_functions/src/cpp/IsEqualTo.cpp
@@ -306,6 +306,9 @@ bool
 isEqualTo(Evaluator* eval, const ArrayOf& A, const ArrayOf& B, bool& needToOverload)
 {
     needToOverload = false;
+    if (A.getElementCount() != B.getElementCount()) {
+        return false;
+    }
     switch (A.getDataClass()) {
     case NLS_DOUBLE: {
         if (A.isSparse()) {