merge of calculateValidSteadyStates and findSteadyState

odesca · Sep 17, 2019 · ed91e9f · ed91e9f
1 parent 95c785d
commit ed91e9f
Show file tree

Hide file tree

Showing 3 changed files with 38 additions and 116 deletions.
diff --git a/BaseClasses/@ODESCA_System/calculateValidSteadyStates.m b/BaseClasses/@ODESCA_System/calculateValidSteadyStates.m
diff --git a/BaseClasses/@ODESCA_System/findSteadyState.m b/BaseClasses/@ODESCA_System/findSteadyState.m
@@ -14,7 +14,7 @@
 %     =====================================================================
 %     name            |  value
 %     ----------------|----------------------------------------------------
-%     method          | 
+%     method          | 'simulate', 'analytically'
 %     inputs          |
 %     simulationTime  | 
 %     showSimulation  |
@@ -119,8 +119,8 @@
                 switch(lower(option))
                     case 'method'
                         value = lower(value);
-                        if ( ~(strcmp(value,'simulate') || strcmp(value,'###ANDRERMETHODE###')) ) %TODO
-                            error('ODESCA_System:findSteadyState:invalidMethod','Method has to be ''simulateFunction'', ''#1'' or ''#2''. The default method was selected.');
+                        if ( ~(strcmp(value,'simulate') || strcmp(value,'analytically')) ) 
+                            error('ODESCA_System:findSteadyState:invalidMethod','Method has to be ''simulate'' or ''analytically''. The default method was selected.');
                         else
                             method = value;
                         end
@@ -208,6 +208,38 @@
         % Return the final value of x
         x0 = val_x(end,:)';
 
+    case 'analytically'    
+        if ~isempty(sys.u)
+            if u0
+                u_steady = u0;
+            else
+                u_steady = sym('u_s',[length(sys.u),1]);
+            end
+            eqns = subs(sys.f,sys.u,u_steady) == 0;
+        else
+            eqns = sys.f == 0;
+        end
+
+        if ~isempty(sys.p)
+            % get parameter
+            temp = sys.getParam();
+            paramVal = sym('p',size(temp));
+            for num = 1:numel(temp)
+                paramVal(num) = temp{num};
+            end
+            % set parameter
+            eqns = subs(eqns,sys.p,paramVal);
+        end
+
+        sys.validSteadyStates = solve(eqns,sys.x,'ReturnConditions',true);
+        x0 = [];
+        for i=1:length(sys.x)
+            eval(['x0 = [x0; sys.validSteadyStates.x',num2str(i),'];']);
+        end
+        if u0
+            x0 = double(x0);
+        end
+
     otherwise
         warning('No correct method selected!');
         x0 = [];

diff --git a/Examples/Systems/PipeSystem.m b/Examples/Systems/PipeSystem.m
@@ -51,16 +51,7 @@
 % Input values for steady state: u0 = [Temperatur In, Massflow In]
 u0 = [40; 0.1];
 % Solve the system equations for the states at the given input values
-steadystate = vpasolve(subs(PipeSys.calculateNumericEquations,PipeSys.u,u0),PipeSys.x);
-x0(1,1) = steadystate.x1;
-x0(2,1) = steadystate.x2;
-x0(3,1) = steadystate.x3;
-x0(4,1) = steadystate.x4;
-x0(5,1) = steadystate.x5;
-x0(6,1) = steadystate.x6;
-x0(7,1) = steadystate.x7;
-x0 = double(x0);
-
+x0 = PipeSys.findSteadyState('method','analytically','inputs',u0);
 ss1 = PipeSys.createSteadyState(x0,u0,'ss1');
 
 %% Linear approximation
@@ -83,27 +74,11 @@
 % Create more steady states and plot a bode plot with all steady
 % states of the system:
 u0_2 = [40; 0.2];
-steadystate = vpasolve(subs(PipeSys.calculateNumericEquations,PipeSys.u,u0_2),PipeSys.x);
-x0_2(1,1) = steadystate.x1;
-x0_2(2,1) = steadystate.x2;
-x0_2(3,1) = steadystate.x3;
-x0_2(4,1) = steadystate.x4;
-x0_2(5,1) = steadystate.x5;
-x0_2(6,1) = steadystate.x6;
-x0_2(7,1) = steadystate.x7;
-x0_2 = double(x0_2);
+x0_2 = PipeSys.findSteadyState('method','analytically','inputs',u0_2);
 ss2 = PipeSys.createSteadyState(x0_2,u0_2,'ss2');
 
 u0_3 = [40; 0.25];
-steadystate = vpasolve(subs(PipeSys.calculateNumericEquations,PipeSys.u,u0_3),PipeSys.x);
-x0_3(1,1) = steadystate.x1;
-x0_3(2,1) = steadystate.x2;
-x0_3(3,1) = steadystate.x3;
-x0_3(4,1) = steadystate.x4;
-x0_3(5,1) = steadystate.x5;
-x0_3(6,1) = steadystate.x6;
-x0_3(7,1) = steadystate.x7;
-x0_3 = double(x0_3);
+x0_3 = PipeSys.findSteadyState('method','analytically','inputs',u0_3);
 ss3 = PipeSys.createSteadyState(x0_3,u0_3,'ss3');
 
 % Linearize all steady states and create a bodeplot