My Blog

CMOD5模型的matlab实现

2022-11-01

%% 数据准备
clear
clc
%% 模型系数准备
global c
c=zeros(1,28);
c(1:13)=[-0.318004,-1.416264,0.006174,-0.058233,0.001851,0.023727,0.071528,0.016858,4.889050,-0.541012,-1.762694,0.487292,-0.058692];
c(14:18)=[0.128558,0.010419,0.08296,0.032500,16.730787];
c(19:28)=[1.466281,7.358231,13.328972,-9.012029,0.999416,2.52457,-0.524396,0.407859,2.227098,-0.306554];
v=0:0.1:20;
theta=30;psai=0;p=1.6;
sigma0=zeros(2,length(v));
for i=1:length(v)
    sigma0(1,i)=sigma0_hh(v(i),psai,theta,p);
end
c(1:13)=[-0.688,-0.793,0.338,-0.173,0,0.004,0.111,0.0162,6.34,2.57,-2.18,0.4,-0.6];
c(14:18)=[0.045,0.007,0.33,0.012,22];
c(19:28)=[1.95,3,8.39,-3.44,1.36,5.35,1.99,0.29,3.8,1.53];
%网站上Fortran程序的参数
%cmod5n(5.0,0.0,45.0) =   0.00872497 这是程序给出的一个参考结果
%c=[-0.6878, -0.7957,  0.3380, -0.1728, 0.0000,  0.0040, 0.1103, 0.0159,  6.7329,  2.7713, -2.2885, 0.4971, -0.7250, 0.0450, 0.0066,  0.3222,  0.0120, 22.7000, 2.0813,  3.0000, 8.3659, -3.3428,  1.3236,  6.2437,  2.3893, 0.3249,  4.1590, 1.6930];
for i=1:length(v)
    sigma0(2,i)=sigma0_vv(v(i),psai,theta,p);
end
res=sigma0_vv(5,0,45,p);
%% 修饰区
plot(v,sigma0(1,:),'b',v,sigma0(2,:),'r','LineWidth',1.5);
xlabel('wind speed(m/s)');
ylabel('\sigma_0');
legend('CMOD5-HH','CMOD5-VV')
figure
plot(v,std2dB(sigma0(1,:)),'b',v,std2dB(sigma0(2,:)),'r','LineWidth',1.5);
xlabel('wind speed(m/s)');
ylabel('\sigma_0(dB)');
legend('CMOD5-HH','CMOD5-VV')
figure
plot(v,std2dB(sigma0(2,:)./sigma0(1,:)),'LineWidth',1.5);
xlabel('wind speed(m/s)');
ylabel('PR(dB)');
%% 核心函数
function res=sigma0_vv(v,psai,theta,p)
res=B0(v,theta)*(1+B1(v,theta)*cosd(psai)+B2(v,theta)*cosd(2*psai))^p;
end

function res=sigma0_hh(v,psai,theta,p)
res=(B0(v,theta)*(1+B1(v,theta)*cosd(psai)+B2(v,theta)*cosd(2*psai)))^p;
end

function res=std2dB(sigma0)
res=10*log10(sigma0);
end
%% B0函数群
function res=B0(v,theta)
x=(theta-40)/25;
res=(10^(a0(x)+a1(x)*v))*(f(v,theta)^gamma(x));
end

function res=f(v,theta)
x=(theta-40)/25;
s=a2(x)*v;
s_0=s0(x);
if s<s_0
    res=((s/s_0)^alpha(s_0))*g(s_0);
else
    res=g(s);
end
end

function res=g(s)
res=1/(1+exp(-s));
end

function res=alpha(s0)
res=s0*(1-g(s0));
end

%% 只和x有关的函数
function res=a0(x)
global c
temp1=c(1:4);
temp2=ones(4,1);
for i=2:4
    temp2(i)=temp2(i-1)*x;
end
res=temp1*temp2;
end

function res=a1(x)
global c
res=c(5)+c(6)*x;
end

function res=a2(x)
global c
res=c(7)+c(8)*x;
end

function res=gamma(x)
global c
temp1=c(9:11);temp2=ones(3,1);
for i=2:3
    temp2(i)=temp2(i-1)*x;
end
res=temp1*temp2;
end

function res=s0(x)
global c
res=c(12)+c(13)*x;
end

function res=v0(x)
global c
temp1=c(21:23);temp2=ones(3,1);
for i=2:3
    temp2(i)=temp2(i-1)*x;
end
res=temp1*temp2;
end

function res=d1(x)
global c
temp1=c(24:26);temp2=ones(3,1);
for i=2:3
    temp2(i)=temp2(i-1)*x;
end
res=temp1*temp2;
end

function res=d2(x)
global c
res=c(27)+c(28)*x;
end
%% B1
function res=B1(v,theta)
global c
x=(theta-40)/25;
temp1=c(14)*(1+x)-c(15)*v*(0.5+x-tanh(4*(x+c(16)+c(17)*v)));
temp2=1+exp(0.34*(v-c(18)));
res=temp1/temp2;
end

function res=B2(v,theta)
x=(theta-40)/25;
res=(-d1(x)+d2(x)*v2(v,theta))*exp(-v2(v,theta));
end

function res=v2(v,theta)
global c
x=(theta-40)/25;
y=(v+v0(x))/v0(x);
y0=c(19);n=c(20);
a=y0-(y0-1)/n;
b=1/(n*((y0-1)^(n-1)));
if y<y0
    res=a+b*((y-1)^n);
else
    res=y;
end
end

关于sentinel-2数据处理的若干技巧

写在前面

随着九月份数学建模国赛的结束，我的竞赛生涯也算是告一段落了。
虽然最后获得了省二的成绩，但是我已经不打算参加了。
一是因为这个比赛过于纯粹，二是也确实没时间了。往后我的精力会主要放在科研上。

1、数据处理的目的

由于我主要是做微波遥感这块，所以对于sentinel-2的处理会较为粗糙。
我主要是想通过sentinel-2图像寻找加拿大北部海区的融池，并且将sentinel-2的数据导出为nc文件供后续处理。
数据处理使用SNAP 9.0。

2、数据读取

可以采用两种方式：直接打开zip文件和打开xml文件。

2022-05-26

rmarkdown测试

前言

　　Matlab 统计工具箱中有 27 种概率分布，工具箱对每一种分布都提供 5 类函数，其命令的字符是：
pdf 概率密度； cdf 分布函数； inv 分布函数的反函数；stat 均值与方差； rnd 随机数生成
　　这里主要介绍常用的rnd类函数 ## rand 函数

rand(n) %返回一个 n×n 在区间 (0,1) 内均匀分布的随机数矩阵
rand(m,n) %返回一个 m×n 在区间 (0,1) 内均匀分布的随机数矩阵

unifrnd 函数

unifrnd(a,b,m,n) %返回一个 m×n 在区间 (a,b) 内均匀分布的随机数矩阵

randn 函数

randn(n) %返回一个 n×n 标准正态分布的随机数矩阵
randn(m,n) %返回一个 m×n 标准正态分布的随机数矩阵

exprnd 函数

exprnd(lambda,m,n) %返回一个均值为lambda的 m×n 的指数分布随机数矩阵

normrnd 函数

normrnd(mu,sigma,m,n) %返回一个满足N(mu,sigma)的 m×n 的正态分布随机数矩阵

poissrnd 函数

poissrnd(lambda,m,n) %返回一个参数为lambda的 m×n 的泊松分布随机数矩阵

binornd 函数

binornd(n,p,m,n) %返回一个满足B(n,p)的 m×n 的二项分布随机数矩阵

randperm 函数

randperm(n) %返回行向量，其中包含从 1 到 n 没有重复元素的整数随机排列
randperm(n,k) %返回行向量，其中包含在 1 到 n 之间随机选择的 k 个唯一整数(1-n随机选择k个互不相同的整数)

perms 函数

perms(n) %返回的矩阵包含了向量 v 中元素按字典序反序的所有排列。  
% P 的每一行包含 v 中 n 个元素的一个不同排列。矩阵 P 具有与 v 相同的数据类型，包含 n! 行和 n 列。

附录

　　以下介绍几种绘制数据分布图的方法，建议配合以上的随机数生成函数学习：
### hist 函数

[N,X] = hist(Y,M); 
%得到数组（行、列均可）Y 的频数表。
%它将区间[min(Y),max(Y)]等分为 M 份（缺省时M 设定为 10）。
%N 返回 M 个小区间的频数，X 返回 M 个小区间的中点。

cdfplot 函数

cdfplot(X) %作样本X（向量）的经验累积分布(也就是F(X))函数图形
%如果表头(title)出现矩形则是中文显示异常，可以通过设置字体解决
%set(0,'defaultAxesFontName', 'Microsft YaHei UI');实测有效，就是每次启动都要设置

参考文献以及资料

1、司守奎数学建模算法课件（2版含源程序）

2、https://max.book118.com/html/2017/0702/119736142.shtm

3、南信大的数学建模课件（Monte-Carlo 方法介绍及其建模应用）

2022-05-26

function varargout = fig2(varargin)
% FIG2 MATLAB code for fig2.fig
%      FIG2, by itself, creates a new FIG2 or raises the existing
%      singleton*.
%
%      H = FIG2 returns the handle to a new FIG2 or the handle to
%      the existing singleton*.
%
%      FIG2('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in FIG2.M with the given input arguments.
%
%      FIG2('Property','Value',...) creates a new FIG2 or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before fig2_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to fig2_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help fig2

% Last Modified by GUIDE v2.5 20-May-2022 19:49:20

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @fig2_OpeningFcn, ...
                   'gui_OutputFcn',  @fig2_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% --- Executes just before fig2 is made visible.
function fig2_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to fig2 (see VARARGIN)

% Choose default command line output for fig2
handles.output = hObject;
global DS1102U
DS1102U = visa('NI', 'USB0::0x1AB1::0x0588::DS1EV193800455::INSTR'); 
fopen(DS1102U); %打开已创建的 VISA 对象 [此处为相关功能语句] 
% Update handles structure
guidata(hObject, handles);

% UIWAIT makes fig2 wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = fig2_OutputFcn(hObject, eventdata, handles) 
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;


% --- Executes on button press in pb1.
function pb1_Callback(hObject, eventdata, handles)
% hObject    handle to pb1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global DS1102U
fprintf(DS1102U,':RUN');
fprintf(DS1102U,'*IDN?');%设置命令 
IDN=fscanf(DS1102U);
set(handles.ed1,'String',IDN);

% --- Executes on button press in pb2.
function pb2_Callback(hObject, eventdata, handles)
% hObject    handle to pb2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global DS1102U
fprintf(DS1102U,':STOP'); 


% --- Executes on button press in pb3.
function pb3_Callback(hObject, eventdata, handles)
% hObject    handle to pb3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global DS1102U
fprintf(DS1102U,':AUTO');
N=DS1102U.InputBufferSize;
fprintf(DS1102U,':wav:data?');
[data,len]=fread(DS1102U,N); 
wave=data(12:end); 
plot(handles.axes1,wave);
% 获取 VPP 
fprintf(DS1102U,':MEASure:VPP?'); 
vpp=fscanf(DS1102U);
set(handles.ed3,'String',vpp/2);
% 获取采样频率 
fprintf(DS1102U,':MEASure:FREQuency?'); 
Fs=fscanf(DS1102U);
set(handles.ed2,'String',Fs);

function ed1_Callback(hObject, eventdata, handles)
% hObject    handle to ed1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of ed1 as text
%        str2double(get(hObject,'String')) returns contents of ed1 as a double


% --- Executes during object creation, after setting all properties.
function ed1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ed1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end



function ed2_Callback(hObject, eventdata, handles)
% hObject    handle to ed2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of ed2 as text
%        str2double(get(hObject,'String')) returns contents of ed2 as a double


% --- Executes during object creation, after setting all properties.
function ed2_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ed2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end



function ed3_Callback(hObject, eventdata, handles)
% hObject    handle to ed3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of ed3 as text
%        str2double(get(hObject,'String')) returns contents of ed3 as a double


% --- Executes during object creation, after setting all properties.
function ed3_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ed3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


% --- Executes on button press in pb4.
function pb4_Callback(hObject, eventdata, handles)
% hObject    handle to pb4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global DS1102U
fclose(DS1102U);
delete(DS1102U);


% --- Executes during object creation, after setting all properties.
function axes1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to axes1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: place code in OpeningFcn to populate axes1

契机

在做导弹拦截这题的时候发现自己不太会用lower_bound和upper_bound,所以写一篇文章备忘。

常规用法

1
2
3

//以下都需要在不下降序列中用
lower_bound(f,f+n,a[i]); //返回f中大于等于a[i]的第一个数
upper_bound(f,f+n,a[i]); //返回f中大于a[i]的第一个数

在不上升序列中用

1
2
3

lower_bound(f,f+n,a[i],greater<int>()); //返回f中大于等于a[i]的第一个数
upper_bound(f,f+n,a[i]，greater<int>()); //返回f中大于a[i]的第一个数
//当然greater也可以换成自己写的比较函数

前言

之前一直考虑复盘一下今年的蓝桥杯，但是苦于没有时间（摸鱼摸掉了）。所以拖到现在才写。
我参加的是C++ B组，最后省一中等偏上。
由于我语文不太好，行文整体偏向于流水账。

Day0

敲了一些常用的数据结构，比如线段树、树状数组、ST表、倍增LCA等等。还写了一下图论最短路的模板题（dijkstra、SPFA、Floyd）。
晚上吃完饭突然发现不能用无线鼠标，于是去买了一个有线鼠标，手感一般（和机房的类似）。买完鼠标后晚上写了一些数论的板子（线性筛、exgcd、线性求逆元等等）。最后看了一些洛谷上一些经典的线性dp（但是没有看递推题单，结果G题推错方程当场炸掉20分）。看完发现也十二点多了，遂睡觉。
晚上比较紧张睡不着，三点钟看了一次手机，大概四点才睡的。

Day1 上午

早上起的比较晚，于是没吃早饭直接开打了。
开考后没几分钟就收到了5张黄牌，吓得我都不敢在动了。但是好在监考老师说没事。于是开始写题。以下是考试的时间线（不准确）：
9:05 写掉A
9:06 看到B和日期有关，遂不写，开始看C
9:14 发现C是模拟，于是切掉
9:15-9:20 开始看D，把玩了一些小案例。发现就是对于每个元素在2(n-i)和2(i-1)之间取最大值。
9:31 写掉D，开始看E
9:40 E看不懂，开始看F，一眼看出暴力前缀和可以得到70。于是直接开始写。
9:50 码出了F的二维前缀和，放着开始看G
9:55 发现G好像有点像原题，通过分析小样例可以一眼看出边界条件，于是开始重点突破递推方程。
10:10 感觉递推方程的f[i][j]意义不明，遂放着，开H
10:15 看完H感觉是搜索，但是又怕爆复杂度，想着重新优化一下降到O(nlogn)以下。开始犯困了。

matlab第四次作业

2022-05-04

第一题代码

clc
clear
a=input('输入三角形的一条边 ');
b=input('输入三角形的一条边 ');
c=input('输入三角形的一条边 ');
if(a>0&&b>0&&c>0)
    if((a+b>c)&&(a+c>b)&&(b+c>a))
        p=(a+b+c)/2;
        s=sqrt(p*(p-a)*(p-b)*(p-c));
        fprintf('三角形面积为%f\n',s);
    else
        fprintf('不能构成一个三角形\n');
    end
else
    fprintf('不能构成一个三角形\n');
end

第二题代码

clc
clear
n=input('输入一个数n ');
sum=0;
for ii=1:n
    sum=sum+ii;
end
sum

第三题代码

function res=myfactorial(n)
res=1;
if n==1
    res=1;
else
    res=myfactorial(n-1)*n;
end
end

第四题代码

a=input('输入一个数a ');
b=input('输入一个数b ');
t=a;
a=b;
b=t;
fprintf('a=%d\n',a);
fprintf('b=%d\n',b);

第五题代码

clc
clear
a=input('输入二次项系数a ');
b=input('输入一次项系数b ');
c=input('输入常数c ');
r=roots([a b c]);
disp(r);

第六题代码

function [r,theta]=Coord_Trans(x,y)
r=sqrt(x*x+y*y);
theta=atan(y/x);
end

第七题代码

clc
clear
n=input('学生成绩: ');
switch floor(n/10)
    case {0,1,2,3,4,5}
        disp('不及格');
    case {6,7}
        disp('及格');
    case {8}
        disp('良好');
    case {9,10}
        disp('优秀');
    otherwise
        disp(-1);
end

第八题代码

clc
clear
a=1;b=1;c=2;n=3;
while c<=500
    c=a+b;
    %fprintf('%d %d\n',a,b);
    a=b;
    b=c;
    n=n+1;
end
fprintf("%d %d\n",c,n);

第九题代码

clc
clear
n=[1 2 3 5 7];
buf=n;
for ii=1:5
    cnt=0;
    while(1)
        if (cnt~=0&&n(ii)==1)
            break;
        else
            if mod(n(ii),2)==0
                n(ii)=n(ii)/2;
            else
                n(ii)=n(ii)*3+1;
            end
            cnt=cnt+1;
        end
    end
    fprintf('%d经过%d步变为1\n',buf(ii),cnt);
end

2022-04-23

第二章

2022-04-23

第一题

　　只需要把$x_1x_2$换成y在利用0-1规划的性质对条件进行变形即可。

第二题

　　很明显决策变量是小区，要达到全部覆盖只需要满足每个小区能选的小学数量大于1即可。
　　这题连数据都不需要。

model:
sets:
col/1..6/:x;
endsets
min=@sum(col:x);
x(1)+x(2)+x(3)>1;
x(2)+x(4)>1;
x(3)+x(5)>1;
x(4)+x(6)>1;
x(5)+x(6)>1;
x(1)>1;
x(2)+x(4)+x(6)>1;
end

第三题

　　没啥好说的引入0-1的限制就行。

model:
sets:
she/1..6/;
com/1..4/;
link(she,com):c,x;
endsets
data:
c=4 2 3 4
6 4 5 5
7 6 7 6
7 8 8 6
7 9 8 6
7 10 8 6;
enddata
max=@sum(link:c*x);
@for(she(i):@sum(com(j):x(i,j))=1);
@for(com(j):@sum(she(i):x(i,j))>1);
@for(link:@bin(x));
end

第四题

　　待更新

第五题

　　首先考虑到比较小的数据范围，可以搜索出料头小于1米的切割方案（大于1米则可以切出一根1米的棍子）。
　　对于这七种切割方案进行规划，即可得到最小使用根数（注意材料根数要是整数）。

model:
sets:
row/1..3/;
col/1..7/:x;
link(row,col):a;
endsets
data:
a=1 2 0 0 0 0 1
0 0 3 2 1 0 1
4 1 0 2 4 6 1;
enddata
min=@sum(col:x);
@for(row(i):@sum(col(j):a(i,j)*x(j))>100);
@for(col:@gin(x));
end

第六题

　　和第三题很像

model:
sets:
car/1..8/;
pos/1..5/;
link(pos,car):c,x;
endsets
data:
c=30 25 18 32 27 19 22 26
29 31 19 18 21 20 30 19
28 29 30 19 19 22 23 26
29 30 19 24 25 19 18 21
21 20 18 17 16 14 16 18;
enddata
min=@sum(link:c*x);
@for(pos(i):@sum(car(j):x(i,j))=1);
@for(car(j):@sum(pos(i):x(i,j))<1);
@for(link:@bin(x));
end

第一章

2022-04-23

第一题

　　没啥好说的，直接用linprog干就完事。只要注意把maxz取负求最小值即可。大于等于的条件也要记得加负号。

clc
clear
c=[3 -1 -1]';
a=[1 -2 1;4 -1 -2];
b=[11;-3];
aeq=[-2 0 1];beq=1;
[x,fval]=linprog(-c,a,b,aeq,beq,zeros(3,1));
fval=-fval

model:
sets:
col/1..3/:c,x;
row/1..2/:b;
link(row,col):a;
endsets
data:
c=3 -1 1;
a=1 -2 1 4 -1 -2;
b=11 -3;
enddata
max=@sum(col:c*x);
@for(row(i):@sum(col(j):a(i,j)*x(j))<b(i));
-2*x(1)+x(3)=1;
end

第二题

　　只要注意把$\vert x \vert$变换为$u_i+v_i$，把$x$变换为$u_i-v_i$进行建模即可。

clc
clear
c=1:4;
c=[c,c]';
aeq=[1 -1 -1 1;1 -2 1 -3;1 -1 -2 3];
beq=[0;1;-0.5];
aeq=[aeq,-aeq];
[uv,fval]=linprog(c,[],[],aeq,beq,zeros(8,1));
x=uv(1:4)-uv(5:8)
fval

model:
sets:
col/1..4/:c,x;
row/1..3/:b;
links(row,col):a;
endsets
data:
c=1 2 3 4;
a=1 -1 -1 1 1 -1 1 -3 1 -1 -2 3;
b=0 1 -0.5;
enddata
min=@sum(col:c*@abs(x));
@for(row(i):@sum(col(j):a(i,j)*x(j))=b(i));
@for(col:@free(x)); !x的值可正可负
end

第三题

　　这里我是用的非线性的优化函数fmincon，但是用线性优化函数再加上整数限制算可能更好。

clc
clear
f=@(x) -(x(1)+x(2)+1.65*x(8)+2.3*x(9)-(300/6000)*(5*x(1)+10*x(6))-(321/10000)*(7*x(2)+9*x(7)+12*x(9))-(250/4000)*(6*x(3)+8*x(8))-(783/7000)*(4*x(4)+11*x(9))-(200/4000)*x(5));
a=zeros(5,9);
a(1,1)=5;a(1,6)=10;
a(2,2)=7;a(2,7)=9;a(2,9)=12;
a(3,3)=6;a(3,8)=8;
a(4,4)=4;a(4,9)=11;
a(5,5)=7;
b=[6000;10000;4000;7000;4000];
aeq=[1 1 -1 -1 -1 0 0 0 0;0 0 0 0 0 1 1 -1 0];
beq=[0;0];
[x,fval]=fmincon(f,randn(9,1),a,b,aeq,beq,zeros(9,1));
fval=-fval

第四题

　　本题使用matlab求解很不方便，这里只提供lingo代码。
　　这道题关键在于对三个仓载货量和最大容许量成比例的建模。

model:
sets:
wu/1..4/:a,b,c,y;
cang/1..3/:w,v;
links(wu,cang):x;
endsets
data:
c=3100 3800 3500 2850;
b=480 650 580 390;
a=18 15 23 12;
w=10 16 8;
v=6800 8700 5300;
enddata
max=@sum(wu(i):c(i)*@sum(cang(j):x(i,j)));
@for(wu(i):@sum(cang(j):x(i,j))<a(i));
@for(cang(j):@sum(wu(i):x(i,j))<w(j));
@for(cang(j):@sum(wu(i):b(i)*x(i,j))<v(j));
@for(cang(j):@sum(wu(i):x(i,j))/w(j)=(@sum(wu(i):x(i,1))/w(1)));
@for(wu(i):y(i)=@sum(cang(j):x(i,j)));
end

!a(i)为运输货物i的重量
!b(i)为单位货物所占的空间
!c(i)为单位货物的利润
!y(i)为四种物资的量
!w,v分别是空间和体积限制

第五题

　　这题和上一题一样，由于决策变量是二维变量使用lingo求解更为方便

model:
sets:
row/1..3/;
col/1..4/:b;
link(row,col):c,x;
endsets
data:
b=62 83 39 91;
c=132 0 100 103
  0 91 100 100
  106 89 100 98;
@text('ex.txt')=@table(x); !存到文本文件中;
enddata
min=@sum(link:c*x);
@for(col(j):@sum(row(i):x(i,j))=b(j));
end

关于sentinel-2数据处理的若干技巧

写在前面

1、数据处理的目的

2、数据读取

rmarkdown测试

cpy

2022/5/26

前言

unifrnd 函数

randn 函数

exprnd 函数

normrnd 函数

poissrnd 函数

binornd 函数

randperm 函数

perms 函数

附录

cdfplot 函数

参考文献以及资料

契机

常规用法

在不上升序列中用

前言

Day0

Day1 上午

第一题代码

第二题代码

第三题代码

第四题代码

第五题代码

第六题代码

第七题代码

第八题代码

第九题代码

第一题

第二题

第三题

第四题

第五题

第六题

第一题

第二题

第三题

第四题

第五题