SNL-WaterPower / WecOptTool-MATLAB

WEC Design Optimization Toolbox
GNU General Public License v3.0
12 stars 9 forks source link

WaveBot Case C #163

Closed ssolson closed 4 years ago

ssolson commented 4 years ago

Description

Example of Multiobjective Function WaveBot Case C converted to use the new file structure.

Fixes # (issue)

Checklist:

ssolson commented 4 years ago

@ryancoe I think the 3D graph looks better with a contoured mesh and using linear axes. image

ssolson commented 4 years ago

In 2D I like labeled contours, with scatter points colored by the contours with the y-scale (Volume) as log. I also reversed the direction of the colorbar so that the color matched the direction of the plot.

I think the contour lines help show the large region where zmax is approximately the same (for large volumes) and further clarifies that along higher value constant power lines as the volume decreases the device must rapidly increase zmax.

image

ryancoe commented 4 years ago

@ssolson can you post the .fig files for me to take a look at? Maybe a .mat file of the results also?

ssolson commented 4 years ago

Here you go. Do you plan to run the code later? plots.zip

ryancoe commented 4 years ago

Any ideas what's going wrong here:

>> ver('MATLAB')
-----------------------------------------------------------------------------------------------------
MATLAB Version: 9.8.0.1417392 (R2020a) Update 4
MATLAB License Number: 10848
Operating System: Mac OS X  Version: 10.14.6 Build: 18G6020 
Java Version: Java 1.8.0_202-b08 with Oracle Corporation Java HotSpot(TM) 64-Bit Server VM mixed mode
-----------------------------------------------------------------------------------------------------
MATLAB                                                Version 9.8         (R2020a)
>> !git show --oneline -s
WARNING: terminal is not fully functional

-  (press RETURN)

ab7a151 (HEAD -> botCaseC, ssolson/botCaseC) Remove comment (cleanup)

>> WaveBot_caseC
Error using simulateDevice>pseudoSpectralControl
Too many input arguments.

Error in simulateDevice (line 13)
            performance = pseudoSpectralControl(motion, varargin{:});

Error in WaveBot_caseC>myWaveBotObjFun (line 182)
    simRes = simulateDevice(deviceHydro,               ...

Error in WaveBot_caseC>@(x)myWaveBotObjFun(x,w,SS,zmax,folder.path) (line 91)
[x,fval,exitflag,output,residuals] = paretosearch(@(x) myWaveBotObjFun(x,w,SS,
zmax,folder.path),nvars,...

Error in globaloptim.internal.FcnEvaluator (line 155)
                fval = feval(FUN, X0);

Error in globaloptim.paretosearch.coldstart (line 50)
    globaloptim.internal.FcnEvaluator(objfun,nvars,[],[],[],[],[],[],nonlcon,optimState.X0,options,"rowwise");

Error in globaloptim.paretosearch.initialize (line 19)
        globaloptim.paretosearch.coldstart(objfun,linConstr,nonlcon,optimState,options);

Error in globaloptim.paretosearch.driver (line 28)
    globaloptim.paretosearch.initialize(objfun,linConstr,nonlcon,optimState,options);

Error in paretosearch (line 231)
[X,FVAL,EXITFLAG,OUTPUT,CINEQ,CEQ] =
globaloptim.paretosearch.driver(objfun,nonlcon,optimState,linConstr,options,OUTPUT);

Error in WaveBot_caseC (line 91)
[x,fval,exitflag,output,residuals] = paretosearch(@(x) myWaveBotObjFun(x,w,SS,
zmax,folder.path),nvars,...

Caused by:
    Failure in user-supplied objective evaluation. Cannot continue.

>> 
ryancoe commented 4 years ago

Also, can you upload the .mat file of the results?

ssolson commented 4 years ago

Also, can you upload the .mat file of the results?

Its in the zip file

ssolson commented 4 years ago

Nope, it all works for me.

>> ver('MATLAB')
-----------------------------------------------------------------------------------------------------
MATLAB Version: 9.8.0.1417392 (R2020a) Update 4
MATLAB License Number: 707564
Operating System: Microsoft Windows 10 Home Version 10.0 (Build 19041)
Java Version: Java 1.8.0_202-b08 with Oracle Corporation Java HotSpot(TM) 64-Bit Server VM mixed mode
-----------------------------------------------------------------------------------------------------
MATLAB                                                Version 9.8         (R2020a)
>> !git show --oneline -s
ab7a151 Remove comment (cleanup)

image

ssolson commented 4 years ago

@ryancoe to follow up on this everything ran to completion without error. I have also attached the final workspace. image botCaseC.zip

ryancoe commented 4 years ago

I am not a huge fan of the contour and surface because of how the interpolation/extrapolation can make it deceptive as to where the space is valid (one may assume that some location on the contour is reachable, but that isn't always true).

What about something like this (marker size is proportional to Zmax)? Would like to add a size legend, seems a bit tricky (https://www.mathworks.com/matlabcentral/answers/228667-how-do-i-create-a-legend-that-explaines-scatter-plot-by-marker-size), but can be done...

image
figure
hold on
scatter(fval(:,1),fval(:,2),fval(:,3) ./ max(fval(:,3))*500,...
    'filled','MarkerEdgeColor','k','MarkerFaceAlpha',0.25)
grid on
xlabel('Power')
ylabel('Vol fun')
ssolson commented 4 years ago

I think that it should be intuitive that the space is not continuous and any point not currently evaluated would need to be evaluated to determine the validity and Obj Fun value of the solution at that point. I think what we have done is a sweep of the space that would allow the designer to narrow in on a region of interest in the Pareto front for their device/ resource and rerun the optimization within that smaller region of interest. I dislike the different size points, while technically quantitative the human abilitity to distinguish between different sizes makes it more qualitative. The same could be said for a continuous contour, which could be discretized into smaller bins (e.g. 10 color bins instead of 256), which is why I added the contour lines. @H0R5E makes fantastic graphics for the papers I've been involved with. Any thoughts on representing the space?

ryancoe commented 4 years ago

hmmm... I'd make the argument that absolute values in this plot are rather meaningless, so it's not too important that people cannot assign a number to a size.

Regardless, I am happy to have @H0R5E take a crack at this.

H0R5E commented 4 years ago

I found a paper that discussed making such plots, and for anything 3D or below, they seem to suggest that scatter plots are the best bet. I would avoid colouring, as it may bias one variable over another. I also think it's important to have the values, because you can identify impracticable results (like the results that tend to zero in this case). I don't think it's an easy task, and I have no great insight, TBH. The best I could come up with is plotting from two difference viewpoints, thus:

ax1 = subplot(1,2,1);
scatter3(ax1, vol, zmax, -pBar, 'filled');
ax2 = subplot(1,2,2);
scatter3(ax2, vol, zmax, -pBar, 'filled');
set(ax2, 'Ydir', 'reverse')
set(ax2, 'Xdir', 'reverse')
set(ax2, 'Projection','perspective');
set(ax1, 'Projection','perspective');
xlabel(ax2, 'WEC Volume [m$^3$]', 'interpreter','latex')
ylabel(ax2, 'PTO Stroke [m]', 'interpreter','latex')
zlabel(ax2, 'Useful Power [W]', 'interpreter','latex')
xlabel(ax1, 'WEC Volume [m$^3$]', 'interpreter','latex')
ylabel(ax1, 'PTO Stroke [m]', 'interpreter','latex')
zlabel(ax1, 'Useful Power [W]', 'interpreter','latex')

pareto

ryancoe commented 4 years ago

Will do something like Fig 3 (without the contour lines) and Fig 4

91176421-7bd58200-e69f-11ea-9c69-100f803e75c7

ryancoe commented 4 years ago

Ok, here's what I came up with, I'm going to cry if you don't like it ;)

WaveBot_caseC_results

fig = figure();
fig.Position = fig.Position .* [1 1 1.5 1]*1;

tiledlayout(3,4,'TileSpacing','compact','Padding','compact')
axb = nexttile([3,2]);

hold on
grid on
scatter(axb, pBar,vol,75,zmax,...
    'filled',...
    'MarkerEdgeColor','k',...
    'MarkerFaceAlpha',0.5);

xlabel('Avg. power, $ \bar{P}$ [W]', 'interpreter','latex')
ylabel('Vol. fun, $(r_0 + r)^3$ [m$^3$]', 'interpreter','latex')

cb = colorbar;
cb.Label.Interpreter = 'latex';
cb.Label.String = ('Max. PTO stroke, $z^{\textrm{max}}$ [m]');
cb.Location = 'northoutside';
set(cb,'YDir','reverse')

set(gca,'Yscale','log')

ylim([min(vol), Inf])

% Create textarrow
annotation(fig,'textarrow',[0.0813492063492065 0.113095238095238],...
    [0.81825396825397 0.554761904761905],'TextEdgeColor',[0 0 0],...
    'TextBackgroundColor',[1 1 1],...
    'String',{'Smaller stroke,','larger vol.'},...
    'HorizontalAlignment','center');

% Create textarrow
annotation(fig,'textarrow',[0.0726190476190479 0.0821428571428572],...
    [0.276984126984128 0.104761904761905],'TextEdgeColor',[0 0 0],...
    'TextBackgroundColor',[1 1 1],...
    'String',{'Larger stroke,','smaller vol.'},...
    'HorizontalAlignment','center');

rlbs = {'$\bar{P}$','$(r_0 + r)^3$','$z^{\textrm{max}}$'};
xlbs = {'Outer radius, $r$ [m]','Max. PTO force, $F_u^{max}$ [N]'};

for kk = 1:2*3
    ax(kk) = nexttile();
end
ax = reshape(ax,[2,3]);

for ii = 1:size(x,2)
    for jj = 1:size(fval,2)
        hold on
        scatter(ax(ii,jj),x(:,ii),fval(:,jj),[],zmax,...
            'filled',...
            'MarkerEdgeColor','k',...
            'MarkerFaceAlpha',0.25);

        if jj ~= size(fval,2)
            set(ax(ii,jj),'XTickLabel',[])
        end

        if ii == 1
            ylabel(ax(ii,jj),rlbs{jj}, 'interpreter','latex')
        end
    end
    xlabel(ax(ii,jj),xlbs{ii}, 'interpreter','latex')
end