Example Trade Study with MATLAB

Information

 
Notes/Code
function [finalResults] = tradstudyMatlab(ConstellationDesignClass)
%tradstudyMatlab Attaches to an open instance of STK10 and does a
%simple constellation design varying the numer of planes and number of
%satellites per plane and return the metrics of choice. A table and 3D
%graph of the results will also generated
%
%NOTE: This will compute coverage based on your parallel computing settings
%in Edit > Preference > Parallel Computing. Here you can force coverage to
%always compute in parallel
%
%   results = tradstudyMatlab(ConstellationDesignClass) takes a class of properties that
%   describes the parmeters you want to vary and the results you desire and
%   runs the analysis in an open STK 10 scenario. This function returns a
%   matrix of the results. If there is an invalid walker configuration then
%   the results will return -9999 for that run
%
%     Example 1:
%     First create a empty ConstellationDesignClass class and then populate the properties  
%
%     parameters = ConstellationDesignClass;
%     parameters.seedPath = '*/Satellite/Satellite1';
%     parameters.deltaWalker = true;
%     parameters.satsRange = 4:6;
%     parameters.planeRange = 2:4;
%     parameters.raanSpread = 360;
%     parameters.innerPlaneSpacing = 1;
%     parameters.coverageName = 'CoverageDefinition1';
%     parameters.useSatelliteAsset = true;
%     parameters.fomPath = '*/CoverageDefinition/CoverageDefinition1/FigureOfMerit/FigureOfMerit1';
%     parameters.dataProviderName = 'Overall Value';
%     parameters.dataProviderElementArray = {'Minimum'; 'Maximum'; 'Average'}; 
%     results = tradstudyMatlab(parameters);

%------------------------------------------------------------------------------
    %Attached to a running STK10 scenario
    uiapp = actxGetRunningServer('STK10.application');
    root = uiapp.Personality2;
    scen = root.CurrentScenario;
    %Convert units to EpSec to be able to compare time values
    root.UnitPreferences.Item('DateFormat').SetCurrentUnit('EpSec');
    %Get path to seed satellite
    seed = root.GetObjectFromPath(ConstellationDesignClass.seedPath); 
    seedName = seed.InstanceName;
    %Get path to coverage definition
    cov = root.GetObjectFromPath(ConstellationDesignClass.coverageName); 
    %Get path to FOM
    fom = root.GetObjectFromPath(ConstellationDesignClass.fomPath); 
    %Set the type of Walker constellation
    if ConstellationDesignClass.deltaWalker
       type = 'Delta';
    else
        type = 'ModDelta';
    end
    
    %Calculate the number of iterations 
    planeCount = length(ConstellationDesignClass.planeRange);
    satsCount = length(ConstellationDesignClass.satsRange);
    raanCount = length(ConstellationDesignClass.raanSpread);
    innerPlaneCount = length(ConstellationDesignClass.innerPlaneSpacing);
    
    %varyParameters and properties are used to find out what the design
    %parameters are and later used for graphing
    varyParameters = [];
    if planeCount > 1
        varyParameters = [varyParameters 1];
    end
    if satsCount > 1
        varyParameters = [varyParameters 2];
    end
    if raanCount > 1
        varyParameters = [varyParameters 3];
    end
    if innerPlaneCount > 1
        varyParameters = [varyParameters 4];
    end
    properties = {'planeRange','satsRange','raanSpread','innerPlaneSpacing'};
    %Number of iterations
    loopNumber = planeCount * satsCount * raanCount * innerPlaneCount;
    loopCounter = 1;
    %Output Array
    finalResults = [];
    invalidRuns = [];
    
    %Show progress bar
    h = waitbar(0,'Calculating...');
    %Nested loop to vary parameters, create walker and query results
    for p = ConstellationDesignClass.planeRange
        for s = ConstellationDesignClass.satsRange
            for r = ConstellationDesignClass.raanSpread
                for spacing = ConstellationDesignClass.innerPlaneSpacing
                    results = [p s r spacing];
                    if spacing >= p
                        invalidRuns = [invalidRuns loopCounter];
                        for i = 1:length(ConstellationDesignClass.dataProviderElementArray)
                            value = -9999;
                            results = [results value];
                        end
                    else
                        %Add Walker Satellites
                        %<NumPlanes> <NumSatsPerPlane> {<InterPlaneSpacing> | <TrueAnomaly>} <RAANSpread>
                        root.BeginUpdate;
                        root.ExecuteCommand(['Walker ' seed.Path ' ' type ' ' num2str(p) ' ' num2str(s) ' ' num2str(spacing) ' ' num2str(r) ' Yes SetUniqueNames']);
                        root.EndUpdate;
                        %Get all satellite in scenario
                        sats = scen.Children.GetElements('eSatellite');
                        satCounter = sats.Count;
                        %Analyzer Coverage Results
                        cov.ClearAccesses;
                        cov.AssetList.RemoveAll;
                        for i = 0:(satCounter-1), % loop through objects in sats
                            currSat = sats.Item(int32(i));
                            flag = strncmpi(seedName, currSat.InstanceName, length(seedName));
                            if (flag && ~strcmp(seedName, currSat.InstanceName))
                                %Add the first child of each satellite as the asset
                                if ConstellationDesignClass.useSatelliteAsset
                                    cov.AssetList.Add(currSat.Path);
                                else
                                    cov.AssetList.Add(currSat.Children.Item(int32(0)).Path);
                                end
                            end
                        end
                        %Compute coverage and get data provider results
                        cov.ComputeAccesses;
                        fomDP = fom.DataProviders.Item(ConstellationDesignClass.dataProviderName);
                        
                        if strcmp(fomDP.Type, 'eDrFixed')
                            fomDP.AllowUI = true;
                            fomDpResult = fomDP.ExecElements(ConstellationDesignClass.dataProviderElementArray);
                            %Add results to temp array
                            for i = 0:(fomDpResult.DataSets.Count - 1)
                                value = cell2mat(fomDpResult.DataSets.Item(i).GetValues);
                                if length(value) == 1
                                    results = [results value];
                                else
                                    results = [results mean(value)];
                                end
                            end
                        elseif strcmp(fomDP.Type , 'eDrTimeVar')
                            fomDP.AllowUI = true;
                            fomDpResult = fomDP.ExecElements(cov.Interval.Start, cov.Interval.Stop, 60, ConstellationDesignClass.dataProviderElementArray);
                            %Add results to temp array
                            for i = 0:(fomDpResult.DataSets.Count - 1)
                                %Get the average value in the array for
                                %output
                                value = mean(cell2mat(fomDpResult.DataSets.Item(i).GetValues));
                                results = [results value];
                            end
                        elseif strcmp(fomDP.Type, 'eDrIntvl')
                            disp('Interval Data Providers not yet supported.');
                            close(h);
                            return
                        else
                            disp('Data Provider not yet supported.');
                            close(h);
                            return
                        end
                        
                        %Unload Walker Satellites unless its the last run
                        root.BeginUpdate;
                        if loopCounter ~= loopNumber
                            for j = 0:(satCounter-1), % loop through objects in sats
                                currSat = sats.Item(int32(j));
                                flag = strncmpi(seedName, currSat.InstanceName, length(seedName));
                                if (flag && ~strcmp(seedName, currSat.InstanceName))
                                    currSat.Unload
                                end
                            end
                        end
                        root.EndUpdate;
                        %Final Output
                        finalResults = [finalResults; results];
                        waitbar(loopCounter/loopNumber, h, ['Run ' num2str(loopCounter) '/' num2str(loopNumber) ' Complete.']); 
                        loopCounter = loopCounter + 1;  
                    end     
                end
            end
        end
    end
    close(h);
    
    f = figure('Position', [100 100 710 410]);
    colnames = {'NumPlanes' 'Sats/Plane' 'RAANSpread' 'InnerPlaneSpacing'};
    columnformat = {'numeric', 'numeric', 'numeric', 'numeric'};
    for i = 1:length(ConstellationDesignClass.dataProviderElementArray)
        colnames{end+1} = ConstellationDesignClass.dataProviderElementArray{i,1};
        columnformat{end+1} = 'numeric';
    end
    
    uitable(f,  'Data', finalResults, ...    
                'ColumnName', colnames', ...
                'ColumnFormat', columnformat,...  
                'Position', [10 10 700 400]);
            
    %Create Title for graphs
    fomType = fom.DefinitionType(4:length(fom.DefinitionType));
    if isprop(fom.Definition,'ComputeType')
       computeType = fom.Definition.ComputeType;
       titleGraph = [fomType  '-' computeType(2:length(computeType))]; 
    else
        titleGraph = fomType;
    end
        
    %If 1 parameter was varied than create a line graph
    if length(varyParameters) == 1
        for i = 5:length(colnames)
            f = figure;
            plot(finalResults(:,varyParameters(1)) ,finalResults(:,i));
            title(titleGraph);
            set(gca,'XTick', ConstellationDesignClass.(properties{varyParameters(1)}));
            xlabel(colnames{varyParameters(1)});
            ylabel([ConstellationDesignClass.dataProviderName ' - ' colnames{i}]);
        end
    elseif length(varyParameters) == 2
        %Needs to be at least 3 points to generate a trimesh graph
        try
            finalResults(invalidRuns,:) = [];
            if loopCounter > 2
                %For each metric create a trimesh and configure
                for i = 5:length(colnames)
                    f2 = figure('Colormap', jet);
                    X = finalResults(:, varyParameters(1));
                    Y = finalResults(:, varyParameters(2));
                    Z = finalResults(:,i);
                    tri = delaunay(double(X), double(Y));
                    mesh = trimesh(tri, X, Y, Z,'LineWidth', 2, 'FaceColor', 'interp');
                    colorbar;
                    title(titleGraph);
                    set(gca,'XTick', ConstellationDesignClass.(properties{varyParameters(1)}));
                    set(gca,'YTick', ConstellationDesignClass.(properties{varyParameters(2)}));
                    xlabel(colnames{varyParameters(1)});
                    ylabel(colnames{varyParameters(2)});
                    zlabel([ConstellationDesignClass.dataProviderName ' - ' colnames{i}]);
                end
            end
        catch
           disp('Not enough unique points to generate graph.'); 
        end 
    end             
end